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Handbook Series Linear Algebra 


Singular Value Decomposition and Least Squares Solutions* 
Contributed by 


G. H. GoLuB** and C. REINSCH 


1. Theoretical Background 
1.1. Introduction 


Let A be a real m Xm matrix with m =n. It is well known (cf. [4]) that 


A=" (1) 
where 
Oey Ve re ST end 2 diaeG a) 


n 


The matrix U consists of ” orthonormalized eigenvectors associated with the 
n largest eigenvalues of A A’, and the matrix V consists of the orthonormalized 
eigenvectors of A’A. The diagonal elements of Y are the non-negative square 
roots of the eigenvalues of A? A; they are called singular values. We shall assume 


that 
0220, 2 oO, = 


Thus if rank(A) =r, 0,.,=0,.2= +-- =o,=0. The decomposition (4) is called 
the singular value decomposition (SVD). 


There are alternative representations to that given by (1). We may write 


tee (=) V? with U?U,=I, 


or 


A= 27 awh ZU eet, and 2) diag (o,5::.40,). 


We use the form (1), however, since it is most useful for computational purposes. 


If the matrix U is not needed, it would appear that one could apply the usual 
diagonalization algorithms to the symmetric matrix A’ A which has to be formed 
explicitly. However, as in the case of linear least squares problems, the com- 


* Editor’s note. In this fascicle, prepublication of algorithms from the Linear 
Algebra series of the Handbook for Automatic Computation is continued. Algorithms 
are published in ALGor 60 reference language as approved by the IFIP. Contributions 
in this series should be styled after the most recently published ones. 

** The work of this author was in part supported by the National Science Founda- 
tion and Office of Naval Research. 
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putation of A? A involves unnecessary numerical inaccuracy. For example, let 


4g 
AZ NEmOe 
0 B 
then 
: 14+p? 1 
a 
so that 


0, (A) = (2+ ?)#, o,(A) = |B]. 


If p?< , the machine precision, the computed A’ A has the form j , and 
the best one may obtain from diagonalization is 6,=|/2, é,=0. ie 

To compute the singular value decomposition of a given matrix A, Forsythe 
and Henrici [2], Hestenes [8], and Kogbetliantz [9] proposed methods based 
on plane rotations. Kublanovskaya [10] suggested a QR-type method. The 
program described below first uses Householder transformations to reduce A to 
bidiagonal form, and then the QF algorithm to find the singular values of the 
bidiagonal matrix. The two phases properly combined produce the singular value 
decomposition of A. 


1.2. Reduction to Bidiagonal Form 
It was shown in [6] how to construct two finite sequences of Householder 
transformations 
PO = T— 2x) 9? (k= 1, 2,..., 2) 
and 
QM —J—2yy? (k=4, 2,...;m—2) 


(where x(*)? x(#) — y(*)? y(*) —4) such that 


PH)... PO AQM ,., Q-9) = EO aes 
@) _ . ey 
Qn 
0 \ im —n) xn 


an upper bidiagonal matrix. If we let A® = A and define 
A+) — Pl) 4) (Reso) Os, 
AG+)) — ACHNQM (k= 14,2, ...,n—2) 
then P) is determined such that 
AG TT =O. (tae tM, ck, 108) 


and Q™ such that , ; 
Gye) =0) (j= h.- 2aeeene 


Singular Value Decomposition and Least Squares Solutions. G. H. Golub etal. 405 


The singular values of J) are the same as those of A. Thus, if the singular 
value decomposition of J®=GEHt 


AS PCy ATO! 
so that U= PG, V=Q4H with P= P®... P™, Q=QM .... Q"—), 


then 


1.3. Singular Value Decomposition of the Bidiagonal Matrix 
By a variant of the QR algorithm, the matrix J is iteratively diagonalized 


so that FEES) pc a Geere » 


where JED = SOTTO TH, 


and S®, T are orthogonal. The matrices T) are chosen so that the sequence 
M) = J* 7 converges to a diagonal matrix while the matrices S are chosen 
so that all J“) are of the bidiagonal form. In [7], another technique for deriving 
{S) and {T} is given but this is equivalent to the method described below. 


For notational convenience, we drop the suffix and use the notation 
J=]®, Taj, s=s0, T=T°, M=jTJ, M=J"J. 


The transition ]-+/ is achieved by application of Givens rotations to J altern- 
ately from the right and the left. Thus 


GS me hmes ji acces rh ml bo Fee bs (2) 
St r 
where 
(e—1)(&) 5 
2 
1 
cos 8, —sin 6, (k —1) 
= sin@, cos 6, (R) 
1 
0 eas 
Oy 


and 7, is defined analogously to S, with , instead of 6,. 
Let the first angle, g,, be arbitrary while all the other angles are chosen so 
that J has the same form as J. Thus, 


T, annihilates nothing, generates an entry {/},,, 
Sf annihilates {J},,, generates an entry {/},3, 
T, annihilates {J},;, generates an entry {J}so, (3) 


and finally 
(See Fig. 1.) 


28* 


ST annihilates {/},,,-1, and generates nothing. 
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— aan 
nie ic” 


Fig.1 


This process is frequently described as “‘chasing”’. Since J =S7 JT, 
M=J'J=T'MT 


and M is a tri-diagonal matrix just as M is. We show that the first angle, gp, 
which is still undetermined, can be chosen so that the transition M—WM is a 
QR transformation with a given shift s. 


The usual QR algorithm with shifts is described as follows: 
M—sI=T,R, 
GF (4) 
T,+sl= M, 


where T,7T,=I and R, is an upper triangular matrix. Thus M,=T,’ MT,. It 
has been shown by Francis [5] that it is not necessary to compute (4) explicitly 
but it is possible to perform the shift implicitly. Let ZT be for the moment an 
arbitrary matrix such that 


toa lex (4) 2 ae), 
(i.e., the elements of the first column of 7, are equal to the first column of T) and 
TTT=I1. 
Then we have the following theorem (Francis): If 
i) M=T'MT, 
ii) M is a tri-diagonal matrix, 
ili) the sub-diagonal elements of M are non-zero, 


it follows that M=DM,D where D is a diagonal matrix whose diagonal elements 
are +1. 

Thus choosing 7, in (3) such that its first column is proportional to that of 
M —s1I, the same is true for the first column of the product T = doyle insihh a Were ns 
therefore is identical to that of 7,. Hence, if the sub-diagonal of M does not 
contain any non-zero entry the conditions of the Francis theorem are fulfilled 
and T is therefore identical to 7, (up to a scaling of column +1). Thus the 
transition (2) is equivalent to the Q R transformation of J’ J with a given shift s. 

The shift parameter s is determined by an eigenvalue of the lower 2 x2 minor 
of M. Wilkinson [13] has shown that for this choice of s, the method converges 
globally and almost always cubically. 
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1.4. Test for Convergence 

If |e, <6, a prescribed tolerance, then |q,| is accepted as a singular value, 
and the order of the matrix is dropped by one. If, however, | e,| <6 for k-n, 
the matrix breaks into two, and the singular values of each block may be com- 
puted independently. 

If ¢,=0, then at least one singular value must be equal to zero. In the absence 
of roundoff error, the matrix will break if a shift of zero is performed. Now, 
suppose at some stage | Qn | = 5 


At this stage an extra sequence of Givens rotations is applied from the left 
to J involving rows (k, k +1), (k, k +2), ..., (Rk, 2) so that 
€xi1={J}i,241 is annihilated, but {J}, ,+2, (J}e41,, are generated, 
{J }s,4+2 is annihilated, but {J}, 2138, {J}a+2,, are generated, 


dfinally | 
and linally {Thy : is annihilated, and {Dunn is generated ‘ 


The matrix obtained thusly has the form 
(F) 


Gy es 


i Ai 
J= q, 0 (R). 


Onia Tria Fate 


Note that by orthogonality 
Ge + Ohya + +H GSE. 
Thus choosing 6=||J||,, 9 (e), the machine precision) ensures that all 6, are 


less in magnitude than ¢9|J||,,. Elements of J not greater than this are neg- 
lected. Hence / breaks up into two parts which may be treated independently. 


2. Applicability 
There are a large number of applications of the singular value decomposition ; 
an extensive list is given in [7]. Some of these are as follows: 


2.1. Pseudoinverse (Procedure SVD) 
Let A be a real m Xn matrix. An »Xm matrix X is said to be the pseudo- 
inverse of A if X satisfies the following four properties: 
i) AXA=A, 
nH) AAA =A; 
ilk MUA eS AX 
i Ay a XAG 


1V 
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The unique solution is denoted by A*. It is easy to verify that if A=UX Ve. 
then At=V2*U* where 2* = diag (oj) and 

1. \tlog for eo,2-0 
— 6) for 0; = O - 
Thus the pseudoinverse may easily be computed from the output provided by 
the procedure SVD. 


2.2. Solution of Homogeneous Equations (Procedure SVD or Procedure Minfit) 
Let A be a matrix of rank 7, and suppose we wish to solve 


Ax = 0 Mor “t= 71h 


where @ denotes the null vector. 


Let 
WES UT, Uy tec | wand: eV a0, 05, ap ee 


Then since Av;—o6,4, (¢=1, 2, ... %), 


Av,=0 for ¢=72+-1,...,% 
and 4,==0,. 
Here the procedure SVD or the procedure Minfit with 6=0 may be used 
for determining the solution. If the rank of A is known, then a modification of 
the algorithm of Businger and Golub [1] may be used. 


2.3. Solutions of Minimal Length (Procedure Minfit) 
Let 5, be a given vector. Suppose we wish to determine a vector x so that 


|}, — A x= min (5) 


If the rank of A is less than m then there is no unique solution. Thus we require 
amongst all x which satisfy (5) that 
||, min 
and this solution is unique. It is easy to verify that 
ee ADD es ee | ee 
The procedure Minfit with p> 0 will yield V, 2, c,..., c,. Thus the user is able 


to determine which singular values are to be declared as zero. 


2.4. A Generalization of the Least Squares Problem (Procedure SVD) 


Let A be a real m Xm matrix of rank 1 and let b be a given vector. We wish 
to construct a vector x such that 


(A +4A)x=b+Ab 
and 
trace(AA7AA) + K?Ab? Ab=min. 
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Here K>0 is a given weight and the standard problem is obtained for K—>0. 
Introducing the augmented matrices A= (A, Kb) and AA =(4A, KAD) and the 


vector 
ft x 
C= : 
—1/K 


we have to minimize trace(4A’ AA) under the constraint (A++AA)*%=0. For 
fixed ¥ the minimum is attained for dd = — AX x*/x" xX and it has the value 
x" A’ Ax/x" X. Minimizing with respect to ¥ amounts to the computation of the 


smallest singular value of the matrix A and ¥ is the corresponding column of 


the matrix V in the decomposition (1) with proper normalization [3]. 


3. Formal Parameter List 


3.1. Input to Procedure SVD 


m number of rows of A, m=—n. 

n number of columns of A. 

withu true if U is desired, false otherwise. 

withv true if V is desired, false otherwise. 

eps a constant used in the test for convergence (see Sec- 
tion 5, (iii)); should not be smaller than the machine 
precision ¢ , 1.e., the smallest number for which 
1+ €)>1 in computer arithmetic. 

tol a machine dependent constant which should be set 


ali:m,1:n] 


u[1:m,1:n] 


equal to B/e) where f is the smallest positive number 
representable in the computer, see [11]. 
represents the matrix A to be decomposed. 


Output of procedure SVD. 
a vector holding the singular values of A ; they are non- 
negative but not necessarily ordered in decreasing se- 
quence. 
represents the matrix U with orthonormalized columns 
(if withu is true, otherwise u is used as a working 
storage). 


v[iin, 1:7] represents the orthogonal matrix V (if withv is true, 
otherwise v is not used). 
3.2. Input to Procedure Minfit 

m number of rows of A. 

n number of columns of A. 

p number of columns of B, =O. 

eps same as for procedure SVD. 

tol same as for procedure SVD. 


ab[1:max(m,n),1:n+)] 


ablt, 7] represents 4, ,, 151m, 15] S7, 
ab[i,n +7] represents 6; ;,1S71Sm, 1 <j. 
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Output of procedure Minfuit. 


ab(1:max(m,n),1:n+/] ab{t,7] represents ¥,, 1575”, 18 qn, 
ab(t,n-+ 4] represents c; ;,, 151 max(m, n), 1595p 
viz. C=U, B. 

g[1:n] same as for procedure SVD. 


4. ALGOL Programs 
procedure SVD (m,n, withu, withv, eps, tol) data: (a) result: (q, u, v) ; 
value m, n, withu, withv, eps, tol; 
integer m, 1; 
Boolean withu, withv; 
real cps, tol; 
array 4, 9, u,v; 


comment Computation of the singular values and complete orthogonal decom- 
position of a real rectangular matrix A, 


A=Udiag(q)V?, UTU=V'V=I, 


where the arrays a[1:m,1:n], u[1:m,1:n], v[1:2,1:n], g[:n] re- 
present A, U, V,q respectively. The actual parameters corresponding 
to a,u,v may all be identical unless withu = withv = true. In this 
case, the actual parameters corresponding to uw and v must differ. 
m=—n is assumed; 
begin 

integer 7,7, k, 1, 11; 

real ¢, /,2,.4.S.%;V,23 

array ¢(1:7]; 

for 1:=1 step 1 until m do 

for 7:=1 step 1 until 7 do u([7, 7] := a[z, 7]; 


comment Householder’s reduction to bidiagonal form; 


[i — aaa Oe 
for 1:=1 step 1 until 1 do 
begin 


AI ea pe ies — is iy Nad BIE 
for 7:=7 step 1 until m do s:= s+u[j, 7]}2; 
if s< tol then g:= 0 else 
begin 
f:= uf[t,t]; g:= if <0 then sqri(s) else — sqrt (s); 
h:= {xg—s; ult, 7]:= f—g; 
for ;:= 1 step 1 until 1 do 
begin 
See OF 
for k:—17 step 1 until m do s:= s+u[k, 1] xufk, 7]; 
[== s/h; 
for k:=7 step 1 until m do u[k, 7] := uk, 7] +/xu[k, 2] 
end 7 
end s; 
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Gite o: St 0} 
for 7:=/ step 1 until do s:= s+ u[i, j]t2; 
if s< tol then g:= 0 else 
begin 
f:=ult,t+41]; g:= if f<0 then sqri(s) else —sqri(s); 
h:=fxXg—s; wlt,i+1]:=f-g; 
for 7:= 1 step 1 until » do e[7] := wt, j]/h; 
for ;:=/ step 1 until ™ do 
begin 
S:= 0; 
for k:=1 step 1 until n dos:=s+u[j, k] xu[i, k]; 
for k:=/1 step 1 until 1 do w[j, k] := u[j, k] +s xe[k] 
end; 
end s; 
y := abs (q[t]) +-abs(e[t]); if y> x then x:= y 
end; 


comment accumulation of right-hand transformations; 
if withv then for 7 := 7 step —1 until 1 do 
begin 
if g+0 then 
begin 
h:= ult,++4] X g; 
for ;:=/ step 1 until 7 do v[j, 7] := u[z, 7]/h; 
for ;7:=/ step 1 until 1 do 
begin 
So — 05 
for k:=/ step 1 until do s:= s+ u(t, k] xv[k, 7]; 
for k:=/ step 1 until ” do v[k, 7] := v[k, 7] +s xv, 2] 
end j 
end g; 
for 7:=1 step 1 until x do v[7, 7] := v[7, 7] := 0; 
Oe. $) oe A Oo olsen 
end; 


comment accumulation of left-hand transformations; 
if withu then for 1 := n step —1 until 1 do 
begin 
l:=1+1; g:= [2]; 
for 7 :=/ step 1 until x do u[2, 7] := 0; 


if g+0 then 
begin 
he ulist) xe: 


for ;:=/ step 1 until x do 
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begin 
Si OF 
for k:—1 step 1 until m do s:=s+u[k, 7] xu, 7]; 
fi:= sth; 
for k:= 1 step 1 until m do u{[k, 7] := u[k, 7] +f xu[k, 7] 
end 7; 
for 7:=17 step 1 until m do u[j, 7] := u[j, 2]/g 
end g 


else for 7:= 1 step 1 until m do wy, 7] := 0; 
ult,7]:= u[t,7] +1 


endz; 


comment diagonalization of the bidiagonal form; 
eps := eps KX; 
for k:=n step —1 until 1 do 
begin 
test f splitting: 
for /:= k step —1 until 1 do 
begin 
if abs (e[l]) < eps then goto fest { convergence; 
if abs (q{l —1]) S eps then goto cancellation 
end /; 


comment cancellation of e{/] if 1>1; 
cancellation : 
Ce OS al een) 
for 1:=/ step 1 until & do 
begin 
fesasxeliierel el t= eox<elz |: 
if abs (f) < eps then goto test / convergence; 


g:= alt]; hi= qt] := sqrt(fxf+e xg); c:=glh; s:= —f Ih; 
if withu then for 7 := 1 step 1 until m do 
begin 
Wr Ud, di 2 oie igo |; 
ulj,i1]:=yxe+2xs; uj, vcs —yxs4+2xXc 
end 7 
end iz; 


test f convergence: 
zi= q{k]; ifl=k then goto convergence; 


comment shift from bottom 2 x2 minor; 

w= Q(t]; y:= q[k—1]; g:= e[kR—1]; h:= e[k]; 
((y —2) X(y +2) +(g —A) X(g +A))/(2XhXy); g:= sqrt (fxf +4); 
((% —2) X(% +2) +hx(y/(if f< 0 then /—g else f +g) —h))/x; 


| 


I 


jie 
ee 


comment next Q & transformation; 
C= Se= 1; 
for 1:=/-+-1 step 1 until & do 
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begin 
ec reltin Yi== giz); i= s Xo er Oxcas 
e[t —4]:—= 2:=— sqrt ({ xf +Axh); c:= }/z; s:= hz; 
fi:=«*%Xc+exXs; gi= —xXxXS+eXC; hi= yXs; yi= y Xe} 
if withv then for ;:= 1 step 1 until 1 do 
begin 
6 =y0[9;t—4];}. 22=9[7, 2]5 
v(j,¢—41]:= xxe+2xs; v[7,7]:= —xxs+2zXc 
end); 
git — A) a= sqrt({/ xf +h xh); C:= fiz; ss= hie; 
fic Xess KOS es — $s Xe Lexy: 
if wethu then for ;:= 1 step 1 until m do 
begin 
yi=uly,t—1]; 2:= ufy, 2]; 
ulj7,t—1]:= yxc+2xs; ufj,t]:= —yxs+2xc 
end; 
end: 
e[l]:= 0; e[k):=f; g[k]:= x; goto test f splitting; 
CONVEL LENCE « 


if z< 0 then 
begin comment g{f] is made non-negative; 


q[k]:= —2; 
if withv then for ;:= 1 step 1 until 2 do v[j, k] := —v{[j, ’] 
end z 
end k 
end SVD; 


procedure Minfit (m, n, p, eps, tol) trans: (ab) result: (q); 
value m, n, p, eps, tol; 
integer m, 7, p; 
real eps, tol; 
array ab, q; 
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comment Computation of the matrices diag(qg), V, and C such that for given 


real mXn matrix A and mx matrix B 


U; AV = diag (q) and UZ B =C with orthogonal matrices U, and V. 


The singular values and the matrices V and C may be used to de- 
termine X minimizing (1) ||A X — Bll, and (2) |X||p with the solution 


X = V x Pseudo-inverse of diag (q) xC. 


The procedure can also be used to determine the complete solution 


of an underdetermined linear system, i.e., rank(A) =m<n. 


The array g[1:2] represents the matrix diag (g). A and B together 
are to be given as the first m rows of the array ab[1:max(m,n),1:n-+p]. 
V is returned in the first 7 rows and columns of ab while C is returned 


in the last # columns of ab (if > 0); 
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begin 
integer 7, 7, k, 1, /1, n1, np; 
real c, /, 2,1, $, %, V, 2; 
array e(1:7]; 


comment Householder’s reduction to bidiagonal form ; 
gi= 4 i= 0; np = 1 2p; 
for 1:= 1 step 1 until 1 do 
begin 
C= os Soe On es AS 
for 7:=17 step 1 until m do s:= s+ ab[j,7]}2; 
if s< tol then g:= 0 else 
begin 
f:= ab[t,1]; g:= if <0 then sqrt(s) else —sqrt(s) ; 
hi fXe—s; abt; 1)v= f—s:; 
for 7:=/ step 1 until 1p do 
begin 
Se OF 
for k:=7 step 1 until m do s:= s+ ab[k, 1] xab{k, 7]; 
fi sik: 
for k:=1 step 1 until m do ab[k, 7] := ab[k, 7] +f xab{h, 7] 
end 
end s; 
Gites ==) 
if 7 <_m then for 7 :=/ step 1 until 1 do s:= s + ab{i, 7]*2; 
if s< tol then g:= 0 else 
begin 
f:= ab{t,++1]; g:= if f<0 then sqrt(s) else —sqrt(s); 
h:=f{xg—s; abltt,i+41]):=f—-g; 
for 7:=/ step 1 until 7 do e[7] := ab[z, 7]/h; 
for ;7:=/ step 1 until m do 
begin 
§ 322 (8 
for k:= 1 step1 until x dos:= s+ abd{j, k] xab[t, k]; 
for k:=1/ step 1 until ~ do aby, k] := ab{j, k] +s xe[k] 


end 
end s; 
y t= abs(q[t]) + abs(e[t]); if y> x then x:= y 


end 7; 


comment accumulation of right-hand transformations; 
for 1:= 1 step —1 until 1 do 
begin 
if gO then 
begin 
h:= ablt,i+1]xg; 
for 7:=/ step 1 until 7 do ab[j, 7] := ab[i, 7]/h; 
for 7:= 1 step 1 until x do 
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begin 
$i== 03 
for k:=1 step 1 until x do s:= s + ab[i, k] Xab[k, 7]; 
for k:=/ step 1 until do ad[k, 7] := ab[k, 7] +s xab[k, 1] 
end 7 
end g; 
for 7:—/ step 1 until ” do ad[i, 7] := ab[j, 7] := 0; 
Gb [8, 3 te ele] sil tae 
end: 
eps := eps Xx; nl:=n-+1; 
for 7:= m-+1 step 1 until do 
for 7 := nl step 1 until xf do ad{i, 7]:= 0; 


comment diagonalization of the bidiagonal form; 
for k:= n step —1 until 1 do 
begin 
test f splitting: 
for /:= k step —1 until 1 do 
begin 
if abs (e(1]) S eps then goto test f convergence; 
if abs (gi —1]) Seps then goto cancellation 


end /: 


comment cancellation of e[/] 1f J >1; 
cancellation: 
60 Se 
for 1:=/ step 1 until k do 
begin 
fiz s Seley elec lal: 
if abs (f) S eps then goto test f convergence; 


gi gle]: qlt)] t= heme sqrt) xf eX e)s c= glk; s = — fk 
for ;:= nl step 1 until np do 
begin 
yy Sab [LIS 4 | 2 abl 2,4 |; 
ab{ll,7]:=cxy+sxXz; ablt,7]:= —sxXy+cxz 
end 7 
end; 


test f convergence: 
z:= q{k]; ifl=k then goto convergence; 


comment shift from bottom 2 x2 minor; 

:= gf]; y:= q[k—1]; g:= e[k—1]; h:= elk]; 

((y —2) X(y +2) + (g—A) xX (g +h))/(2XAXy); Br sart(fxf+1); 
((x —2) x (x +2) hx (y/(if f<0 then /—g else / +g) —A))/x; 


| 


Ge 
ee 
fe 


I 


comment next QR transformation; 
ci=si=1; 
for 1:=/+1 step 1 until & do 
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begin 
g:= eft]; vyi= alt]; i= sXg; g-— cXE; 
e(t —1]:= 2:= sqrt (fxf+hxh); c:=flz; s:= hz; 


fi=xxXe+eXs; gi —xXKs+ EXC; hi= yXsS; Vi= YX, 
for ;:= 1 step 1 until 7 do 
begin 
x% i= ablj,4—1]; z:= ably, 2]; 
ab[j,t—1]:= *xc+2xs; ablj,1]:= —xxXs+2xC 
end 7; 
q(t —1] := 2:= sqrt({xf+hxh); c:= flz; s:=hlz; 
f:i=cxg+sxXy; x:= —sxXg+exy; 
for ;:= nl step 1 until np do 
begin 
Vy t== abit —1, 9) 5, 2:= able, 7]; 
abli—1,7]:=cxy+sxXz; ablt,7]:= —sxy+ecxzZ 
end 7 
endz; 
e{l]:= 0; e[k]:=f; g[k]:= x; goto test f splitting; 
convergence: 


if z< 0 then 
begin comment q[k] is made non-negative; 


q{k]:= —2; 
for 7:= 1 step 1 until ~ do ab[j, k] := —abdlj, k] 
end z 
end & 
end Minfit; 


5. Organizational and Notational Details 


(i) The matrix U consists of the first » columns of an orthogonal matrix U,. 
The following modification of the procedure SVD would produce U, instead of U: 
After 


comment accumulation of left-hand transformations; 
insert a statement 
if withu then for 1:= +1 step 1 until m do 


begin 
for 7:=n-+1 step 1 until m do u[1, 7]:= 0; 
Wit, 4) == 4 

end 7; 


Moreover, replace 1 by m in the fourth and eighth line after that, i.e., write 
twice for 7 := / step 1 until ™ do. 


(il) m =n is assumed for procedure SVD. This is no restriction; if m<n, 
store A’, ie., use an array at[1:n,1:m] where at[i, j] represents a; ; and call 
SVD (n, m, withv, withu, eps, tol, at, q,v, uv) producing the mxm matrix U and 
the » Xm matrix V. There is no restriction on the values of m and n for the 
procedure Minfut. 
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(iii) In the iterative part of the procedures an element of J“ is considered 
to be negligible and is consequently replaced by zero if it is not larger in magnitude 
than ex where ¢ is the given tolerance and 

x = max (|q;| +|e,|). 


1Sisn 
The largest singular value o, is bounded by x/|/2 So, < x2. 


(iv) A program organization was chosen which allows us to save storage loca- 
tions. To this end the actual parameters corresponding to a and u may be identical. 
In this event the original information stored in a is overwritten by information 
on the reduction. This, in turn, is overwritten by w if the latter is desired. Like- 
wise, the actual parameters corresponding to a and v may agree. Then v is stored 
in the upper part of a if it is desired, otherwise a is not changed. Finally, all 
three parameters a, wu, and v may be identical unless withu = withv = true. 

This special feature, however, increases the number of multiplications needed 
to form U roughly by a factor m/n. 


(v) Shifts are evaluated in a way as to reduce the danger of overflow or under- 
flow of exponents. 


(vi) The singular values as delivered in the array q are not necessarily ordered. 
Any sorting of them should be accompanied by the corresponding sorting of the 
columns of U and JV, and of the rows of C. 


(vii) The formal parameter list may be completed by the addition of a limit 
for the number of iterations to be performed, and by the addition of a failure 
exit to be taken if no convergence is reached after the specified number of itera- 
tions (e.g., 30 per singular value). 


6. Numerical Properties 


The stability of the Householder transformations has been demonstrated by 
Wilkinson [12]. In addition, he has shown that in the absence of roundoff the 
Q Ralgorithm has global convergence and asymptotically is almost always cubically 
convergent. 

The numerical experiments indicate that the average number of complete 
QR iterations on the bidiagonal matrix is usually less than two per singular 
value. Extra consideration must be given to the implicit shift technique which 
fails for a split matrix. The difficulties arise when there are small q,’s or @,’s. 
Using the techniques of Section 1.4, there cannot be numerical instability since 
stable orthogonal transformations are used but under special circumstances there 
may be a slowdown in the rate of convergence. 


7. Test Results 


Tests were carried out on the UNIVAC 1108 Computer of the Andrew R. Jen- 
nings Computing Center of Case Western Reserve University. Floating point 
numbers are represented by a normalized 27 bit mantissa and a 7 bit exponent 
to the radix 2, whence eps = 1.549 —8, tol =4) —31. In the following, computed 
values are marked by a tilde and m(A) denotes max|q; ,|. 
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First example: 


224010 2 3 IZ —1 1 0) 
14 OHO 0 8 2-1 4 
—1 13 —1—11 3 | es Fo ee 2 
ha — Fay AS ee gt Pes 5 et the ; 
) 8 1 2 4 Gh 
9 Anaee Z 5 4 +ah 6 3 
2 —6 6 5 4 tities © oan © 
4 5 0 —2 2 eS. ae 


== V1Q4S, da = 20,-. Og-=VIO4, Og 0,0 


The homogeneous system Ax=6 has two linearly independent solutions. Six 
QR transformations were necessary to drop all off-diagonal elements below the 
internal tolerance 46.4,, —8. Table 1 gives the singular values in the sequence 
as computed by procedures SVD and Minfit. The accuracy of the achieved de- 
composition is characterized by 


~~ 


m(A —USV7) =238,,—8, m(UTU —1 =81,,—8, m(V7V —1) =3.39—8 


Because two singular values are equal to zero, the procedures SVD and Minfit 
may lead to other orderings of the singular values for this matrix when other 
tolerances are used. 


Table 1 
Op O, — 6, 
0,065 7 —9.6 
19.595 916 191 
19.999 999 143 x 1078 
1.9719 —7 —19.7 


35.327 038 518 


The computed solutions of the homogeneous system are given by the first 
and fourth columns of the matrix V (Table 2). 


Table 2 
dy V4 Dy oy Uys, 
— 0.4190 9545 0) Ae 0 (Def.) 
0.4405 0912 0.4185 4806 1 Wy 0.6 
— 0.0520 0457 0.3487 9006 12 —1.3 <107* 
0.6760 5915 0.2441 5305 x ECG) 0.3 


0.4129 7730 — 0.8022 1713 453 — 0.8 
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Procedure Minjfit was used to compute the solutions of the minimization 
problem of Section 2.3 corresponding to the three right-hand sides as given by 
the columns of the matrix B. Table 3 lists the exact solutions and the results 
obtained when the first and fourth values in Table 1 are replaced by zero. 


Table 3 
My Xe #5 % #y 3 
—1/12 ) — 1/12 = 0:0583 3333 0.1749 — 8 — 0.0833 3333 
fe) 0 0 — 0.581) —8 — 1.094) —8 — 1.114) —8 
1/4 ) 1/4 0.2500 0002 1.5519 —8 0.2500 0003 
—1/12 ) —1/12 — 0.0833 3332 0.7419 — 8 — 0.0833 3332 
1/12 (0) 1/12 0.0833 3334 0.3319 — 8 0.0833 3334 
Residual ys 
0 8V5 815 


A second example is the 20 X21 matrix with entries 


(0) au eS3 ; 
A Cady cael St 20 
a; ~= 421 —4 if +=] 
ue Gn | 
—1 if <4 Baie S 


which has orthogonal rows and singular values o,_,=Vk(k +1), k=0, ..., 20. 
Theoretically, the Householder reduction should produce a matrix J with 
diagonal — 20, 0, ..., 0 and super-diagonal —/20, og, ..., Og9. Under the influence 
of rounding errors a totally different matrix results. However, within working 
accuracy its singular values agree with those of the original matrix. Convergence 
is reached after 32 QR transformations and the 6,, k=1,..., 20 are correct 
within several units in the last digit, Gj; =1.61,) — 11. 

A third example is obtained if the diagonal of the foregoing example is 
changed to 

4,;=1, 451520. 


This matrix has a cluster of singular values, 019 to 0;, lying between 1.5 and 1.6, 
O29 = 2, O2,= 0. Clusters, in general, have a tendency to reduce the number of 
required iterations; in this example, 26 iterations were necessary for convergence. 
Gx, =1.49,9 —8 is found in eighteenth position and the corresponding column 
of V differs from the unique solution of the homogeneous system by less than 
3.4,9 —8 in any component. 

A second test was made by Dr. Peter Businger on the CDC 6600. 

A third test was performed on the IBM 360/67 at Stanford University. The 
example used was the 30 x30 matrix with entries 


GC t-i>7 
—{ if 4<}7. 


The computed singular values are given in Table 4. 
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Table 4. Singular values 


18.2029 0555 7529 2200 
2.4904 5062 9660 3570 
1.8059 1912 6612 3070 
1.6253 2089 2877 9290 
1.5546 4889 0109 3720 
1.5217 8003 9063 4950 
1.5060 4262 0723 9700 


6.2231 9652 2604 2340 
2.2032 0757 4479 9280 
1.7411 3576 7747 9500 
1.6018 3335 6666 2670 
1.5440 8471 4076 0510 
1.5166 4741 2836 7840 
1.5038 0424 3812 6520 
0.0000 0000 2793 9677 


3.9134 8020 3335 6160 
2.0191 8365 4054 5860 
1.6923 5654 4395 2610 
1.5828 6958 8713 6990 
1.5352 8356 5544 9020 
1.5123 8547 3899 6950 
1.5021 1297 6754 0060 


2.9767 9450 2557 7960 
1.8943 4154 7685 6890 
1.6547 9302 7369 3370 
1.5673 9214 4480 0070 
1.5279 2951 2160 3040 
1.5088 8015 6801 8850 
1.5009 3071 1977 0610 


1.5002 3143 4775 4370 


Note that 09/0, + 1.53 X1071° so that this matrix is very close to being a 


matrix of rank 29 even though the determinant equals 1. 
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Numerical Integration of Products of Fourier 
and Ordinary Polynomials* 


D. G. BETTIS 


Received June 6, 1969 


Abstract. Sets of coefficients for four finite difference methods of numerical 
integration are presented that will integrate without truncation error products of 
fourier and ordinary polynomials. These sets are formulated such that they are free 
from computational difficulties. 


I. Introduction 


In a previous paper [4] sets of modified integration coefficients for the Cowell 
method of numerical integration of order six were given which had the property 
that they would integrate without truncation error products of a fourier poly- 
nomial and an ordinary polynomial. These sets of coefficients were characterized 
by their explicit formulation. The purpose of this paper is the extension of the 
set of modified Cowell coefficients to any order of the integration method, as 
well as the development of similar sets of modified integration coefficients for 
the St6rmer, Adams-Moulton, and the Adams-Bashforth methods of numerical 
integration. These modified coefficients will be given explicitly in a form such 
that they can be computed from a simple algorithm which does not suffer from 
numerical difficulties. 

The four methods of numerical integration based upon backward differences 
will be referred to as follows [2]: 1) Cowell — the implicit method for differential 
equations of the second order in which the first derivatives are absent; 2) St6rmer 
— the explicit method for these second order differential equations; 3) Adams- 
Moulton — the implicit method for differential equations of the first order; and 
4) Adams-Bashforth — the explicit method for differential equations of the first 
order. 


II. Integration Formulae 


For the integration of a function /(¢), tabulated at equally spaced values of 
the independent variable ¢ with the step-length h, any ascending diagonal of a 
difference table can be used for the development of a method of numerical 
integration [4]. In such a table, let m be an integer, ¢,,=mh, f,,=f(t»); and let 
A° f (m) =n» Af(m =) =n Pe A? (m) =n 1 2hmn i 96 ect. For a dif- 
ferential equation of the form 


PHO) _ #(x, 1) 


* Part of this work was formulated earlier by the author in a dissertation presented 
to Yale University in partial fulfillment for the degree of Doctor of Philosophy. 


29* 
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the double integration may be obtained by 


(1.1) fy —2%y-+ 24 =* Y 4, A*# (p — 5) 
k=0 


where %),=1, P=1 corresponds to the classical Cowell formula, and =O re- 
presents the Stormer formula. 


For the integration of a differential equation of the first order, 


dx (t) 
TO _ 4 (x, 0) 
a suitable integration formula is 
~ h 
(1.2) m= h Ya, AY (63): 


When # =1 the formula is the classical Adams-Moulton method, and when =0 
it is the classical Adams-Bashforth method. Since the sets of «’s are assumed 
to be constant in these classical methods, the solution of the differential equations 
is expressed as a linear combination of the differences of a particular diagonal 
of the difference table. The coefficients « are different for the four methods of 
integration. However, with this distinction well in mind, no ambiguity will result 
if the same symbol is used for all four of the integration methods. 


In order to determine the four sets of coefficients «, assume that 
(2) pa) = 2" 


where Z represents any fixed complex number. Then it follows that for the 


double integration 
(i) a 


. h(1—Z) |? 
jh 
bal drdt=Z| = ie 


0 th 


and for the single integration 


Letting 


the integration formulae (1) become 


6.1) Sar asi 
4 4 ‘ 

(3.2) 1—€ |log(1 = =r \e) 

(3.3) face 


Ue (A loga—H = PC), 
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where P(¢) is defined as 
(4) P(¢) Pte oa Xp oh 


k=0 
and where the integers 1, 2, 3, 4 denote the Cowell, Stérmer, Adams-Moulton, 
and the Adams-Bashforth methods, respectively. 
Upon expanding the left-hand sides of (3) in powers of ¢, with |¢|<1, the 
coefficients of similar terms of € may be compared, yielding, for the first six 
coefficients the following sets [2]: 


Cowell method — 
1 


ti — 1 Kp = 12? =O, 
(5.4) a ant Alert he = 221 
a RAH POE G FOE GG 486 
St6érmer method — 
1 1 
ay = Ox Xs —— 42. ? Xs 42. ? 
ne aee89: be 5 p03 
%4 = 340? = 40? %6 = 42006” 
Adams-Moulton method — 
1 1 1 
pera 2 ee ae as? 
6-3) iat |) Sees, ge 865 
%4= "720? %s = 4440’ %6 = 60 480 ° 
Adams-Bashforth method — 
eee Bie cages 
Fs Fe ont oe = 427 Ga gar? 
(5.4) 255 __ 475 _ 19087 
Xa = 920” ms = 4440.’ %6 = 60480 ° 


The order, 7, of the integration method is defined by the subscript of the 
highest coefficient «, retained in (1). It is important to remember that the classical 
sets of integration coefficients are based upon the expansion of the left-hand 
sides of (3) in an infinite power series in terms of the variable ¢. If, in application, 
(1) are used with the order chosen such that either the highest coefficients are 
zero, or the difference A” of the function being integrated vanishes identically, 
then (1) will integrate the function exactly. However, in practice these two 
conditions are usually not satisfied. In effect, the integration of the function is 
being approximated by the truncated polynomial of the right-hand side of (3), 


(6) Pi(C) =D uC". 

k=0 
Therefore, the integration of a function by the use of (1) with a finite set of 
coefficients will result in an error because the formulae are based upon an ap- 
proximating polynomial instead of an exact polynomial. 
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While formulating a method to reduce these truncation errors for the four 
classical methods of numerical integration when they are applied to functions 
which are products of fourier polynomials and ordinary polynomials, the basic 
formulae (1) will not be altered, but the coefficients « will be modified. 


III. Formulation of the Problem 


Before deriving the modified coefficients which will integrate without trunca- 
tion error the products of fourier and ordinary polynomials, the modified coef- 
ficients will be obtained which will integrate only fourier polynomials exactly. 
The final sets of coefficients which will integrate the products of the fourier and 
the ordinary polynomials will be obtained as limiting cases of these first sets 
of coefficients. 

Much insight into the general problem of deriving the modified coefficients 
results from the determination of the coefficients which will integrate the functions 
sinwt and cos w t. Thus, for (2) let /(¢) =exp (i m ¢), and introducing the notation 
20=wh, Z in (2) becomes 

Li exp(210)y 
and ¢ becomes 

€=1—exp(— 270). 
Therefore, Eqs. (3) reduce to: 


(7) L,(o)=F,(0), 7=1, 2, 3,4, 
where 
sin o \2 
1, (0) = (2°) exp(—2¢0), 


oO 
sin 
L,(c) = ——* exp(—io), 
L, (6) = ~~ exp (io). 


The coefficients « must be chosen such that (7) are satisfied for the given fre- 
quency w, as well as for —qa, if the four integration methods of order are to 
integrate the functions /(¢) =sin w ¢ and f(t) =coswt exactly. 

For the integration methods of order two the coefficients x, and «, must be 
chosen such that the L,(c) are identical to 


1+ 0%¢ + aC, 


or 


1+ a [1 —exp(—270)] +a, [4 —exp(—270)]?. 


For the methods of order two, the modified coefficients can be determined 
by solving the two equations which result from (7) when the two frequencies w 
and —w are considered. If o is not zero, or if, for the two Adams methods, o is 
not an odd integer multiple of z/2, then the solution of the two equations yields 
for the 
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Cowell method — 

(8.1) = —1, > = A(x); 
Stérmer method — 

(8.2) a, = A*(u) —1, Oy = A(x); 
Adams-Moulton method — 

(8.3) % = w*(u) —1, Oy = (ws); 
Adams-Bashforth method — 


(8.4) m= (3—u)wt—1, a =m (u) +u*(u); 
where 
woh ‘ 
eters =: u—4sin2¢, 
and where 
he sin o \2 sin o 
2(u) = ( o ie PMY) acbaee 


A(u) = : (Sars = ay 


w=+(stc-sxt 
ag ~ 4 \sin?o¢ osinocosa)’ 


With these choices of %, and «, the four integration methods of second order wiil 
integrate the functions sin w ¢ and cos w / exactly, except for any round-off error 
which accumulates in the algorithm. 

Henceforth, the coefficients which have been adjusted in the above manner, 
will be referred to as modified coefficients, and the coefficients of (5) will be 
referred to as the classical coefficients. For the modified coefficients two sub- 
scripts, 7 and m, will denote the 7-th coefficient of the set and the order of the 
method, respectively. As o becomes very small, the above coefficients approach 
the corresponding classical coefficients. 

If the order of the integration method is greater than two, then the a; ,, 
({=1, 2,..., m) coefficients can be determined such that the functions cos w ¢ 
and sin wt are integrated exactly. It will be assumed that is even and equal 
to 2v. Thus, there will exist y distinct frequencies @,, @,,...,@, such that the 
2y coefficients will integrate the functions cosm,¢ and sinw,t (k=1, 2,...,”) 
exactly, if (7) are satisfied for the corresponding o, and for —o,. Since the right- 
hand sides of (7) are polynomials in the variable ¢, the above condition is equiv- 
alent to the following problem of polynomial interpolation: 

Given 2y+1 points ¢,, k=—y»,...,—1,0, +1,..., +», where ¢,= 
4 —exp(—27o,), and where the points are on a circle which has a unit radius 
in the complex ¢-plane, construct a polynomial P,(¢), 7 =2y, such that L ;(o,) = 
PHC), 4 =4) 2,3, 44and,where o_, = —o;,. 

Once the coefficients « are determined such that the above conditions are 
satisfied, the integration methods will integrate without truncation error not 
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only the functions sin w,t and cos @,t, but also the linear combination 


(9) Ay + > (a cos w,t +d, sin w,t). 
kel 


The constant a, is integrated exactly because the polynomial P,(¢) begins with 
unity. 


IV. A Special Case 


When only two of the 2 coefficients are modified, the remainder of the coef- 
ficients retaining their classical values, the resulting set of coefficients will inte- 
grate exactly the function 


(10) Q,,-2(t) + acoswt+bsinwt, 


where Q,,(¢) is an ordinary polynomial of degree m in the variable ¢. In this 
special case there are given two points ¢,=1—exp(—270,), k= —1, +1, on 
a unit circle in the complex ¢-plane passing through the orgin, where the remaining 
nm —2 points are located. The problem is the construction of a polynomial P, (¢) 
such that the relations (7) are satisfied for the two values of k. This problem 
can be solved by the technique which was used for the modification of the coef- 
ficients for the integration methods of order two. Thus, the problem reduces 
to the solution of the two equations (k= —1, +1): 


(11) L,(o,) =F, (¢,), += 1, 2, Dacte 


for two of the coefficients, for example «, ,, and « 


n,n nm—1,%° 
Before solving for the two unknown coefficients it is advantageous to intro- 
duce two sequences of polynomials R,,(u) and S,,(u) defined by the recurrence 
relations [4]: 
Rnri= “(Ry — Ry-1), Ry= 2, R,=4, 


(12) 
Smt1= 4(Sm—Sm—-1)> So=0, S,= 4. 


For m=0, 1, 2,..., the following relations can be proved by induction: 


(—41)"F¥sin' 20 So5, 


sin 2mo = 


(13) cl ae 
sin(2m+1)o= Sa es : 

and 

(14) S sin (EA) Pet eet ee OR as 


Rome = (—1)"u"t* (2m +1 + QR), ; 


where QF (w) is a polynomial with constant coefficients of degree # in the variable u 
with the term of degree zero absent. Also the following preliminary relations 
will be needed: 

m diss4 pads T Sin (up) 
cpnepanianiey 


Ch oe = Ry, (u;) , 


, 


(15) 


where t = 27 sin 20. 
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Solving (11) for «, ,, and «,_1, gives 


1 PEE eee Ss ES | — p— (0{-*# —C45* | 
+ 2i%,n (= I =] , 


Lp n= url t 
where %,=1, g=”—1 when p=n, g=n when f=n—41, and L;,=L;(o;). 
Using (15), the definition of ¢,, and the relation, 
exp (1.6,) = (¢,)-1 27 sin o,, 

and defining the first term in the bracket as *;(w), then 

By = A*(u) yea: 

Bx. = —A*(u) S,/u, 

Fs = —m*(u) R 


q-1) 
EX, =*(u) Royal. 
Therefore, if the classical coefficients are used for «,,, J=1,2,...,%—2, 


then the last two coefficients are 


RR eas E* — pad 
Opn = seme — |Fps() + Don? Sy al, 
k=0 


p=n,qaH=n—1; p=n—-1,9¢=nN. 


(16) 


With these coefficients the integration methods will integrate the function (10) 
exactly. 


V. Modified Coefficients for Distinct Frequencies 
To integrate the function (9) without any truncation error requires the solu- 
tion of the problem of interpolation given in Section III. Any of the classical 
interpolation formulae are available for this purpose. Here the technique will 
be to obtain the two highest coefficients by an application of the Lagrangian 
interpolation method: 


: = olin 
(17) L;(o) i IT ee n= 2y. 


The lower coefficients could be obtained by the same technique, but this results 
in a rather tedious procedure that involves computation in the complex variable ¢. 
Therefore, a more efficient algorithm in the real field is developed in which the 
lower coefficients are determined by recurrence with respect to the coefficients 
of lower order. The final goal is the determination of the two highest coefficients 
in explicit form for any order. Thus, after the remaining coefficients are found 
by the recurrence method, the complete set of modified coefficients for the four 
methods will be readily available. 

The derivation of the expression for the highest coefficient for the four 
methods of any given even order will be given. The modified coefficient «, ,, 
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equals the coefficient of the ¢” term of the Lagrangian interpolation formula: 


Lin €—n(C—2=1)" 1 — Ly = Se (Ce —1)”—* ru 
(18) Hells. Up Ms (Cp—G— it f 


eat Dm, 


where 


v 
1 
iho) fF Taree 


m=1 
m=+k 


The expression (18) can be expressed in terms of wu, A(u), and w(u) for the four 
integration methods. 


Since, for the Cowell method, 
U 


Lyn ia re (1—¢,), 
the right-hand side of (18) becomes: 


“1 (1 —fn) C4 (C—p—1)"?— (1 —C_ 2) Se (Ge — 1)? 
a 402 (f,—C_2) ine 


(19) 


Zo 
Since 


jae a gaye Oe 
the numerator of (19) becomes 
C5, (C,—1" CO +oGa—n 
or, because ¢, = 1—exp(27 0,), the numerator reduces to 
(—1)’—" [27 sin 2(y — 2) o,] + (—1)’ [27 sin 2(v —1)o,]. 


Using the definitions (13), the above expression simplifies to 


So(v—2) , S2(v—1) 
| hat + = c 


and, since the denominator is equal to 40?7, (18) becomes 


T* 
403° 


(20) . [=a So(y—1) 
} = 


y—1 
UR Uy 


With the aid of the first formula of (14), the expression in brackets can be shown 
to be of the form 


(—1)’"14+Q% 2(u,)], 


where 


Qm (ut) = Gu + ou? + --- +9,,u™ 


Therefore, by using the finite expansion of 1/(u,u,...u,) from Section VII, (19) 
becomes of the form 


(21) (4) > [Fae see Oyen 6 oe &s 1+9, Up+ st quaauynat n* 


2 
=a) Ur 403 
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Since the constants c are arbitrary, they can be chosen such that the numerators 
of the two terms in the brackets of (24) are identical, i.e., the coefficients of 
the polynomial 


1+ 04% +--+ +0, 0 
are selected so that the polynomial equals 


Se@v—2) , So—1) _ Sav—s 
or a a 
k 


ur Tie 


Thus, (21) can be expressed as 


which is the «, ,, coefficient for the case of distinct frequencies for the Cowell 
method. 


Defining F, ,(u) as 
(22.1) F, (a) = — 252) 10), 


CY hfe 


the «,,, coefficient for the Cowell method becomes the (y—1)-th divided dif- 
ference of the function F, ,(u) at the nodes w,, w.,..., #,. In a similar manner 
it follows that the «,, ,, coefficient for the Stérmer (7 = 2) and the Adams-Moulton 


(t =3) methods is the (vy —1)-th divided difference of the function F, ;(u), where 
(22.2) jag es PH 

Ike l 
(22.3) F, 5(u) = aes) (a), 


Likewise, for the Adams-Bashforth method (a,,,,— 3) is the (vy —1)-th divided 
difference of F, ,(u), where 


(22.4) B, 4(u) = — 72) nu). 


2u’ 


The factor one-half appears in this method because a term of degree (vy —1) is 
necessary in the finite expansion of the first term of (18). 


The second highest coefficient equals the coefficient of the ¢"~* term of the 
Lagrangian interpolation formula: 


a Linon (C—1— 1)? 2 —Ly — abn (Cp 1)"— 
(23) 2 sar Te Up a n (—t k) Uk “IT, 


k=1 ; 
+= 1, 2, 3,4, 


where 
= (1 = » “) Ky 
k=1 


With the aid of (12), (13), and (14), and with a development analogous to that 
for the highest coefficient, («,_ ,, —2*) becomes the (vy —1)-th divided difference 
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of the function F,-1i(¥), tma4, 255, A where 


Boi =— 1M), 
F,42(u) = 2) au), 
2 
i. F,_4,9(u) = A wu), 
F,-1,4() = — a2 p(w). 
In order to obtain the remaining coefficients «,, ,, #7=1, 2,..., 7% —2, assume 


that the set of coefficients of order » —2 are known and that the two highest 
coefficients have been computed. Consider two polynomials P,(¢) and P,_2(¢) 
interpolated at the nodes 


(25) PASTE ST ee ey 
and 
(26) PREC TTA l=vy—1, 


respectively. The polynomial having the points (26) as zeros is 


C (C2 — my + ty) (0? — Ugh + tty) -.. (0? — 04,0 +m). 


This polynomial which vanishes at the » —2 zeros, multiplied by a linear factor, 
is the difference of the two polynomials P,(¢) and P,_ ,(¢), 


(27) F,(6) —F—2(6) =C (0? — myo + my) (C2? — gd + tg) --. (0? — uy, f +m) (ro —S). 


The parameters 7 and s can be obtained as a function of the two highest coef- 
ficients by first expressing the left-hand sides of (27) in terms of the modified 
coefficients of order m and  —2, respectively, and by then comparing the coef- 
ficients of €” and ¢"~", yielding 

l 


Y=On no $= Oy —a,n tT Onn 2 Uz. 
=1 


With these values for 7 and s, (27) becomes 


I 1 
(28) P,(¢) ae 8 (¢) + On,nb a fe Xn—1,0 - n,n 2s Me ut (a ax u,C? ai uC) . 
The lower coefficients can now be obtained by comparing the coefficients of 
similar terms of ¢ in (28). This technique is easily adaptable to a computer sub- 
routine. 


VI. Modified Coefficients for Confluent Frequencies 


The modified coefficients previously developed will integrate without trunca- 
tion error the function (9). However, if two or more of the frequencies w, approach 
a common limit, for example if @,+>@,—--->@,,, the modified coefficients 
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for this case of confluent frequencies will integrate exactly the function 


v 


(29) A +Qn—1(4)| Dd) (a,cos@,t +b, sin w,?)], 


=m-+1 
i.e., the product of an ordinary polynomial and a fourier polynomial. 


In order to obtain the sets of modified coefficients which will integrate the 
above function, the limits must be obtained of the modified coefficients for distinct 
frequencies as the m frequencies approach the value of confluency. Since the 
modified coefficients for distinct frequencies are the divided differences of the 
function F, ;(«), p=," —1, and 1—1, 2, 3,4, at the nodes m,, %,..., 4,, the 
desired limits are available from the well established theory of divided dif- 
ferences [3]. 

Sets of these modified coefficients will be given without proof for three of 
the more important limiting cases: a) two, b) y —1, and c) all of the y frequencies 
approach a common limit. The first case is applicable when the two highest 
coefficients for distinct frequencies suffer from a loss of significant digits because 
the differences of uw, and w, in their denominators become small. As the confluent 
frequencies approach zero in the second case, the coefficients approach those of 
the special case of Section IV. The third case is characterized by the lack of the 
troublesome factors of wu in the denominators of the two highest coefficients. 
Only the two highest coefficients will be given since the lower coefficients may 
be computed by the recurrence technique of Section V. 

In the first case assume that w, and w, approach the common limit w. If 
G,,(#) is defined as 


v 


G, (u) =[] (u—,), n= 6, 


k=3 
then 
Fy; = u 
(30) pn = Ky, + D [eae | + DB, ale) IT p=n,n—I, 


with the definitions 


and 


156 ea” 40 #=al, 2, 3,4, 


n—1,7 


and where D =d/du. For the second case let the first (v —1) frequencies approach 
the common limit w, then 


sis Uy arte Dal bd oes p=n,n—1 
Uy—U)?—1? 


(31) Opn =Ky s+ 7 (u—u = 


where the operator D is defined as 
dk 
duk’ 


Diya Dies ce ee 
For the third case all the frequencies approach the limit w. Here the highest 


coefficients become 


Dia re ) 
(32) bp n= Ky; - Ger S ’ p=N, n— 1. 
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VII. A Finite Expansion 


The development of the two highest coefficients for the case of distinct fre- 
quencies requires an expansion of the function 1/(a#, u,... %,). 


The (vy —1)-th divided difference of the function g(u) =1/u is [3}: 


(33) Fier, teil Stes Ve ate 


where the (vy —1)-th divided difference of a function g(u) is defined by 


(34) fate t= Delos, = TT aay: 


Combining (33) and (34) gives 


(35) = eee pitt 
R=L 


Since the (y—1)-th divided difference of a polynomial of degree (vy —2) is zero, 
then 


ca Uk 


Furthermore, the (vy —1)-th divided difference of a polynomial of degree vy —1 
equals the coefficient of the term of the polynomial of degree » —1, or 


(37) We yt Se. 


fal te 


Adding (35), (36), and (37) yields the desired finite expansion of the function 
1] (Uy Ug... U,) : 


(38) 1g PSN taut toni) tek rm 
ws ; ra a . 
k=1 


where the constants c are arbitrary. 


VIII. Some Power Expansions 


During computation the functions A*(u), A(u), u*(u), and w(u), and their 
derivatives will suffer from a loss of significant digits if the variable o becomes 
small. This loss of precision may be avoided by expanding the functions in a 
power series in terms of the variable w. 


The expansions of A*(w) and A(u) are presented in [4]: 


1 1 2 Bi 


* = = ty 75 Peas 
AN) eal ai sult macea a ta Te hC ONS , 
1 1 yu dly 289 
AWN ieiang Honeeunan 3628800 “ bo 
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In order to obtain similar expansions for u*(w) and «(w) as a function of 
the variable w consider the following integration formula [1]: 


ax (3)=frmaemalr+ 547 (5)— 5 4410) — 3 441(5) 


0 
11 i 1 191 191 1 
+ 420 4°71) + F440 4°7(2) 60480 4°F() — 430 960 Ade r 4]. 


The function to which this formula will be applied is f(#) =cos w ¢. Since t,,=mh, 
and since 2o0=@ h, it follows that /,,=cos 2mc. Also, since 


cos 2(m +1)o —cos2mo= — 2sin(2m+1)osino, 
then 


Af(m) = — 2sin 20 sin2mo — > cos 2mo. 


Furthermore, 


k 
A?*—-1 #(m) = (—1)* |u*®—? sin 20 sin 2m a + “> cos 2mol, 


or, the odd first differences of cosw ft are 


A*-* #(0) = > (—w). 


Likewise, since 
cos 2(m +1)o0 + cos2(m —1)o0=2cos2a0cos2mo, 
then 
A?}(m) = —ucos2mo, 
and, in general, 
A?* #(m) = (—u)*cos2mo. 


Thus, the even central differences of cos wt are 


A®* #(0) = (—w)*. 


From the integral 
h 


[ coswtdt= * sino h, 
0 
it follows that 


1 sin wh h sin o cos o 
Ve nes . 
2, (a) Oo 


With these result (39) becomes 


SEE AS) hod Ce 1 2 1 ; 23 


= a= qt 
o 6 190° Wsi2” 226800 ~ 


Using this relation the expansions of w*(u) and (wu) are: 


1 1 11 191 
* 2 3 ices 
NM) — ae, @ me yag “1 agooeo” 


1 ie, ees Olmee, 2497 og 
1 (%) = — 45 — 520 “ — Go480 3 628 800 
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IX. Elimination of Numerical Difficulties 


The two highest coefficients suffer from a loss of significant digits during 
their computation when the differences of u,, U2, ..., U,, ™ = 2, in the denomina- 
tors become small. This occurs when m of the frequencies approach a common 
limit. This difficulty may be avoided by using the set of coefficients for /+m 
confluent frequencies, where / denotes the number of confluent frequencies of 
the original set of coefficients. For example, if this situation appears because 
@, and w, approach a common value when the modified coefficients for distinct 
frequencies are computed, the set of coefficients (30) should be used. 

When the variable « becomes small the functions A*(u), A(u), u*(u), w(u) 
and their derivatives become indeterminant. If this occurs the power expansions 
from Section VIII will eliminate the difficulty. 


If o=(2m+1)>, 
derivatives, and the derivatives of the functions A*(u) and A(u) become infinite. 
For a given problem the value of w is specified, but the value of # can be varied. 
Therefore, since 20=qwh, a suitable choice of the step-length will avoid the 


problem of o being an odd integer multiple of /2. 


m=0,1,2,..., the functions w*(u) and w(u) and their 
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for the Navier-Stokes Equations 
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Abstract. An application of stability analysis by the energy method is made to 
a practical problem of unsteady viscous flow solved numerically by Fromm. The 
practical stability criterion for the scheme is determined and a rigorous proof of 
convergence to a smooth solution is given. It is also shown how to construct an 
energy for any 3-level scalar difference equation. 


1. Introduction 

The investigation of stability of difference schemes for initial-value problems 
has largely amounted in practice to linearising the perturbation equations about 
a smooth solution, and setting the coefficients constant. The techniques of Fourier 
analysis are then applicable. This heuristic procedure is not always satisfactory 
since the effects of boundary conditions and nonlinearities may become important. 
The energy method can sometimes provide a rigorous analysis of these problems, 
and is complementary to Fourier analysis in the sense that the main difficulty 
in applying it is to construct a suitable energy for the constant coefficient scheme. 
Analysis of variable coefficient and nonlinear schemes may then be comparatively 
straightforward. 

The present paper carries out this program, taking as an example the difference 
schemes used by Fromm [3] to calculate unsteady viscous flows past an obstacle 
in a channel. The construction of an energy for a given difference scheme is an 
unsolved problem in general, but a solution is given for a class of difference 
schemes which includes Fromm’s. This scheme requires the vorticity to be specified 
on the obstacle and walls of the channel. It is shown that the boundary condition 
chosen by Fromm has no adverse effect on stability. Finally, the convergence of 
the finite difference solutions to sufficiently smooth solutions of the Navier-Stokes 
equations is established. 

Most of the work presented in this paper is contained in a D. Phil. thesis 
submitted to the University of Oxford under the supervision of Dr. K.W. Morton, 
to whom I am indebted for many helpful discussions. I am grateful to the Uni- 
versity of Edinburgh and the Science Research Council for providing financial 
support. 

2. Statement of the Problem 

The Navier-Stokes equations for two-dimensional viscous incompressible flow 

may be written in the form: 
Ke) 0 7) 
Bia ee) AaeAs (VQ) =9V22 


PeY+Q=0 


(2.1) 
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where » is the kinematic viscosity, 2 the scalar vorticity, and WY is the stream 
function related to the cartesian components U and V of the fluid velocity by 
U = dW dy, V = — 00x. The specific problem we shall consider is that of un- 
steady flow past a stationary obstacle in an infinite channel OS y<1. The 
boundary conditions are that the relative velocity between the fluid and the 
walls or obstacle is zero. It is assumed that the initial conditions can be adequately 
specified by some irrotational flow and an impulsive translatory motion of walls 
or fluid in the x-direction. We may then expect a wake to develop on the down- 
stream side of the obstacle. 

Fromm [3] has obtained a numerical solution of this problem by finite- 
difference methods. He assumes that the influence of the obstacle on the flow 
is negligible at points sufficiently far downstream, say at x =L. The flow variables 
at points on this line are assumed equal in value to those at corresponding points 
along a line upstream of the obstacle, say at x =0. Thus we impose a periodic 
boundary condition in the x-direction. The main features of the problem are 
shown in Fig. 1. We denote by R the closed region outside and on the obstacle 
I Omilex< | Oni: 


Fig. 4 


It is assumed that the problem has a unique solution such that 2, YeC4(R) 
for 0<¢ST. If this assumption is not true, nonlinear instabilities may develop 
during the computation at points where the solution is varying rapidly, for 
instance round the leading edge of the obstacle. Such effects were indeed observed 
by Fromm in cases where the Reynolds number was high. 


3. The Finite-Difference Equations 

A square mesh with spacing / is set up over R and it is assumed that the 
boundaries of the obstacle pass through the mesh points. Let the mesh point i, 7 
be at (th, 7h) for OSt SI =L/h, OS] SJ =1/h. We define the following opera- 
tors: 

(a) Translations: T.,@;;=@;41,,- 

(b) Differences: 4,,=-+(T,,—J), 

A 


0x — elec) 
(c) Averages: 14, =4 (Ts, +J), 
box = a(R += J) 
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Similar operations are defined in the y-direction. In the difference equations 
we omit subscripts, thus w/; ~Q(ih, jh, n At) is written as w”. w and wp are defined 
at mesh points and w,v as difference quotients of y; for convenience we shall 
choose to define them at mesh points also by 


n - n n = 
u"=htAy,y", v=—h1A,y”. 


All mesh functions are regarded as elements of the real Hilbert space /,(R) 
endowed with inner product 


(¢, 1) = 2:5 hs; and norm ||| =//(¢, 4). 
We also define the discrete Fourier transform 
F ($) = Did; ;exp[—th (hy + 7h). 
R 


Occasionally we shall deal with the space ¢,(R) with norm || =sup | ¢;,|. 
R 


||| and || are related by the inequalities || < ||, |¢| <KA+|¢|| where K 


is a constant depending on R. 


The difference equations used by Fromm to solve (2.1) are the following: 


wo" t* —@"* + 2A(Ao,(u@)” + Ag, (v@)") = 20(u, tuo" —@"**—@""*) (3.44) 


and 
V2y" +o" =0, (3.1b) 


where A=At/h, 0=2y At/h?, and V? is the standard 5-point approximation to 
the Laplacian. In (3.1a) an approximation of the leapfrog type is used for the 
transport terms, and the viscosity terms are approximated by a generalisation 
of the Dufort-Frankel scheme. Eqs. (3.1) are to be solved at all points in the 
interior of R. The boundary conditions are y constant along each wall and on 
the obstacle. We also require the following to hold: 


for allt, %; 9—=U, (velocity of the wall y=0); v;9 = 03 
and u;7= U, (velocity of the wall y=1);  v,; =0; 


also u;; =v, ; =0 for all points on the obstacle. In addition there are the periodicity 
conditions on w and yp: for all 7, @,;;=@o;; Yr; = Yo;- Finally, in order to solve 
(3.1a) we must specify values of w on the boundaries of R; these must be estimated 
by extrapolation, and may produce numerical instabilities. Fromm chose on 
physical grounds to set the boundary values equal to those at a distance h in- 
wards along the normal; we shall show in § 5 that this choice does not affect 
stability. In § 4 we construct an energy for (3.1a) with # and v constant, in § 5 
we show that the practical stability criterion for (3.1a) is A(|u| +]|v|)< 3, 
and in § 6 we prove that solutions of (3.1) converge as O (h?) in f, to a C* solution 
of.@.1) im K. 
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4. Analysis of the Constant-Coefficient Scheme 
We need only consider the one-dimensional scheme 


pipe 


ot? —o" 1 4 26 A," = 0 (24,0 —o” ow") (4.4) 


since the energy for the two-dimensional problem will be an obvious generalisation. 
In [2] an energy for (4.1) was constructed in an intuitive manner from first 
principles. However this is not easy, and we give here a systematic solution of 
the problem which can be applied to any 3-level scalar difference equation. 


We define the vector $”=[w”,w”"~"]’ and take the discrete Fourier trans- 
form ¥ of (4.1). Using the relation 
F (Ts, w) = e*# Fw) 


where £=kh we obtain ¥(¢”"*t1) =GF(o”) where G is the amplification matrix. 


Here 


where «= 2(0cos& —iosiné)/(1+ 0), B=(1—0)/(1+0). 


The stability of (4.1) is defined in terms of G; we require ||G”|| to be uniformly 
bounded in & for all 7. A necessary and sufficient condition for stability which 
we shall use was derived by Kreiss [4]. This is that there must exist a constant 
Cy >0 and for each € a positive definite Hermitean matrix H(&) such that 


CpISHSCyI and G*HG<H. 


(By ASB we mean that A —B is non-positive definite.) 


Now the energy for (4.1) will be a quadratic form in w”, w”~* and powers 
of their translates T,.w”, etc. This corresponds to a matrix H in Fourier space 
whose elements are polynomials in e**. It is an unsolved problem as to whether 
an H of this form exists for every stable scheme. However for matrices 


cf 


a suitable solution is proportional to 


Proof. We have to show that stability >H =>CyiI>0 and G*HG<H. In 
terms of « and f these conditions become 


2(1+|B|?) —|a/?-2C>0 (4.2) 
and 


Best, a tae] <1—[plp (4.3) 
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First we determine the stability conditions on G. Application of the Kreiss- 
Buchanan criteria ([4,4]) on the triangular form of G produces the following 
conditions on the eigenvalues A, and A,: 


(von Neumann condition) 
and 
| Ay| [4+ A, Ag| SM max (| A, — Ag], 
for some M>0. 
Now, since the L.H.S. cannot vanish when there is a double root on the unit 
circle, the latter condition is equivalent to 


1—| Ay |p| |) 


max ( 


ye KS 0. (4.4) 


Using the relations 4,A,=6 and 2|A,|?+2|A,|?=|a|?-+|A,—A,|? we express 
(4.2) as 
| Ax — Ag]? + 2 (4 —|Ay|?) (1 —]Ag|4) 2C >. (4.5) 


To show that (4.4) > (4.5) we note that if the maximum of the L.H.S. of (4.4) 
is |A,—A,|=K then (4.5) is satisfied with C=K*. If on the other hand the 
maximum value is 1—|A,| =< then since |A,—A,|? = (|A,| —|A,|)? the L.H.S. 
of (4.5) is greater than /{(K, uw) =(K —p)?+2K(2—K) w(2—p) where w= 
4 —|A,|. From the fact that f(K, uw) has no stationary value with 0<u<1 for 
fixed K <1, we can show that /(K, 1) =/(K, uw) =f (K, 0) =K*. Similarly for 
the case 1—|A,| = K. Thus in all cases (4.4) implies (4.5) with C = K?. 

Finally, since A.A, =8 we must have |f| <1, and application of the von Neu- 
mann condition leads to 


z|o|?+ 2/6 +202] S1+|6]?. (4.6) 


Isolating |B + ¢?| then squaring and expanding out, we obtain |a+a,|< 
4 —|B|? which is (4.3). 


Remark. H is constructed as follows: 


|B\?=b and |c—fe—apl?S(b —|B|*) (4 —|a|?—2Alac —O). 
b and c are chosen to make the R.H.S. a complete square; thus 
Emimatiy Pmalalbl). 
Using this form for H, the corresponding energy for (4.1) is 
Sn = (1+) (lo"P + lo" P) + 2e|0"P 
—20(1+@) (0, uo @" 7) +20(1+@) (w”, Ao"). 


The difference scheme is stable when S, =C>0 and S,,,<S,. This in turn 
implies conditions on the mesh ratios 9 and o, which can be determined from 
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conditions (4.3). However they can be derived more directly from an analysis 
of the difference equations themselves. This involves some tedious algebra, but 
is essential in any case for the variable coefficient analysis. 


First we define norms of diagonal differences 
Dy= lo" —T, oP + flo" —T_ oh 
and 
F,=4(1 —o) |o” —T, wo 72? +3(1+0) lo” —T_ ow”? fP. 
Writing (4.1) as 
wth tw") = 2u,@" — 3 (w"*? —w"* + 20 A, 0") (4.7) 
and taking the inner product with w”*? —w"—'+204,w” we obtain 
lo” +? |? — jo” "2 + 26(@"*7* +0" 7, Ayo") 
= 2 (wt? —@"*, uy") + 40 (Ug@”, Aw”) (4.8) 
—o7|o"*? —w"* + 264,0"[P. 
Applying (A.6), (A.1) and (A.2) from the appendix we find that 
FE. 


ata — B= —@7|o"** —@"* + 20 Ayo”. 
(Note that F, is not a suitable energy since it is not equivalent to the /,-norm.) 
By taking the inner product of (4.1) with w”t?+ "+ we obtain a second 
identity 
Jon Jo" tp = 200" Lo", ye) 
— 20 (wt? +."-}, Ayo”) — ojo" +o". (4.9) 


From (4.8) and (4.9) we obtain 
oan — De = 20 (w"** si Ce on [y@") a 20 (w"** arr, A,o") 
+20 (0"*, of" * — py") —20(00", ©” — 190") 
— ola"? +o""|P oe o||o"* Sao a i 20 A,w" 2 
+ 20(1+ 0) (o”*?+ 0", 4,0") 


—2007?(w"t* + w""1, Ayo"). 


By separating out the terms independent of o and applying (A.1) we find after 
cancellation and rearrangement that 


Snt1—Sn= — 20 [o"P —2(o"~*, woe”) + lo” P] 
—e(0"*? —@""* + 264)", 20 Ayo") 
—o (or Oa 20 A,w") 


+200(1—0) (o"t1+@", Ao"). 
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Then, on substituting (4.1) in the product on the second line and using (A.1) 
we find that 


Snti—Sp= —2eD, +40Q(0", Ayo”) 
=z —2@ n 


=0 i jo =. 


Finally, to show that S,, is positive definite we apply (A.1) and (A.6) and 
write it in the form 


S,,= lo" +o Ayo"? P + Jo" [2 — 0240" P 
+ @*D, + 0(o"+c4.o"—, wo” —T, wo"? +0" —T_o*). 


Estimation of the inner product by the Cauchy-Schwarz inequality and application 
of (A.4) leads to 


S,, =o" +00" + oD, + (1 — 0) lo" 
—o[e|o" +cA,o"*|? + eD,] 
> (108) fo" Ip, 
Therefore if |o|<1 we have 


lou Cs = C Se ora S,. 


Since S,, is obviously bounded above in terms of the f,-norm, we have shown 
that the scheme is stable in /, if |a|<1. 


5. Variable-Coefficient Analysis and Boundary Conditions 


In this section we derive the practical stability criterion for (3.1a) with u 
and v given and investigate the effects of boundary terms. We have to make 
the a priori assumption that w and v are Lipschitz continuous; this will be justified 
in §6. First we obtain identities analogous to (4.8) and (4.9). (3.1a) can be 
written alternatively as 


wt) ee ae 57% + a,” 


—102[w"*t—@" 1 + 21(A9,(u@)” + Ag, (vo)”)}. (5.1) 
Then, taking the inner product of (5.1) with 
w"tt — o""14 21(Ao,(uw)” + Ao, (v@)") 
leads to 
Ilo” +22 — lo" 2 + 2A(w" +? +0", Ay, (wm)” + Ao, (v@)”) 
= (o"t! _o"", a, +H, 0") rs 


= 2A(Ao,(ua)” alow (vo)", My + Myo") 
— 297 \a"t* — "1 + 24(Ap, (wo)”" + Ao, (v@)”)P. 
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On the other hand, taking the inner product of (3.1a) with w”t*+a"~* yields 
lo" +32 — |o" 12 + 2A(w"*? +04, Ay, (wo)” + Agy (vm)”) 
= 20 (ot? +0}, uw, + yo") (5.3) 
(ae 20|\o"** +o" |P. 
Our object is to show that (3.1a) is stable in /, if max (A| |, A|v|)< }. The proof 
is along the same lines as that for the one-dimensional constant-coefficient scheme, 
but we have to consider extra terms which arise from the variable coefficients, 
and values at the boundaries. The presence of the obstacle is neglected for the 
moment. 
The first type of inner product to arise in (5.2) and (5.3) is 


(w"—*, Ay, (wa)”) = — (w”, u” 4y,.0"—") 
= — (0, u"* Ay, 00""?) + (08, WI — WF Ay 0") 
= — (@", Ap, (u@)"~*) +0 (A) jo” | o” |, 


by (A.7) and the periodicity condition on the vorticity. 
Using (A.7) we also have 


(w"—?, Ao, (vw)”) = — Ox Ao, (v@)"~*) +0 (h) |o"—* | ||o” | 


+¢ Ss hfopzt 1 (vo)? 7 + azz (v@)} pa 


—. ots (vm)i1 —wzz* (va)? o}- 


Since v7 9=vj,7= 0 for all 7, 2 the sum reduces to 


+E (ony! (vm)? 7- 1— Fo ‘(v)i1}. 


Now we shall show in § 6 that |v —V| =O(A) in the interior of R. Since V=0 
on the boundaries and its derivatives are bounded in R, V=O(h) at mesh points 
adjacent to the boundaries and hence v=O(A) also at these points. Thus the 
sum is at most O(h) ||w”~*| lw” |. 


The second type of inner product producing boundary terms is 


(Hoyo ™*, oo”) = (wo, sae oo”) 


+h?/2- ah {ory of 1— OF aa * 07 t7—ori oF oto wr, 

on application of (A.1). If we define w at the boundaries 7=0 and 7=J by 
@? o>=@7, and w} ;=wj 7_1 for all n, we see that the above sum vanishes, so 
that this boundary condition has no adverse effect on stability. This is in fact 
the condition which was used by Fromm. It is evident that the above argument 
remains valid when the obstacle is included. 

The third type of inner product we have to consider is (49, (w@)” + Ap, (vm)", 
Mz + My"). 

By (A.8) this is O(h) |\”|? plus boundary terms in the y-direction. We have 
shown already that these can be neglected. 
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Let us now define the norm of diagonal differences 
n= 40° —T,,0° P+ Zo" —T_,o" 
+4 lo" —Tyo"™P +20" —T_,o" 


=|o"P+[o"P— (0%, a, Fo") 


if w is defined as above at the boundaries, and the energy 


Sn= lo"? + lo" P+ @(@", 20" —p, Fu, 0") 
+ @?D, +2A(1+ @) (@”, do, (u@)""! + Ay, (vw)"-). 
When we calculate S,,,,—S, we perform the same manipulations as in the con- 


stant coefficient analysis and obtain the corresponding result, apart from ad- 
ditional terms which we have shown to be 


20 (h) |o”| lo” | = 0 (4?) (lo" |? + Jo" P) 
when the correct boundary condition is used. We arrive at 
Sn41—S,= —20[D, + 24(0", Ag, (uo)"* + Ag, (v@)"“")] 
+0 (A?) (lo" P+ lo”). 


Using the expanded form of D,, we split up the R.H.S. into terms involving 
translation operators in each dimension separately: Thus 


Se ee) Neca “yy dleee! Oh alan ane-7) EN Os Ee, ot) 
+#2(0* —T_,w", 14+2AT_,u"- wo" —T_,w"—})] 
+40 (h) |o”|| |o"™| 
+ corresponding terms in the y-direction. 
The inner products are positive definite if 2A max(|w|, |v|)<1, and then 
Snat cy Sn = 0 (A t) (lo” P "ts lo" [P). 


Finally we show that S, is equivalent to the /,;-norm, by the same argument 
as for the constant coefficient scheme. Using (A.8), we write 


S,=|lo”* +1449, (uw)""* + AApg, (v@)"*P 
+ lo" — 2740. (woo)”™ + Aoy (vo) P 
+0 (A?) oP + 0?D, 
+ 0(w" + AA, (uw)"* + AAp, (v@)"", 20" — nw, + yo" *) 
=> (1 —O(A2)) |w”“* |? — 224, (uw)”—* + Ay, (v@)”* P 
= (1 —0 (At) —2?(|u| +] 2])?) o”P 
on estimation of the inner product in the same way as before. Hence if 
Amax(|u|, |v|)<4 then a norm equivalent to the f,-norm has a growth which 


depends only on the Lipschitz constants of w and v. Thus the practical stability 
criterion for (3.1a) is A max(||,|v|)< #. 


444 N. G. Campbell: 


6. Convergence to a Smooth Solution 


In this section we prove that the solutions of (3.1) converge as O (h*) in f,(R) 
to a C4 solution of (2.1). First we obtain the error equations by substituting 


a=2 +a’ 
u=U+y 
Vea. 
pat +e 


in (3.1). We have 
wt wy!" 4 21 [Ag, (Ww 2 +Uo' +u'w’)” + Ag, (v'.24+Vo'+0'o’)"] 
=20(u, +u,o'"—o'"t*—o'") + ALT, 
Yay Oo =F, (5.4b) 


(5.4a) 


where 7), and 7), are the truncation errors of (3.1). It can be verified by Taylor 
expansions that if Q, Y are C4 in R, 

| Z,,| =O( (At)? +h? + (At/h)?) 

| Z| =O.(7). 
For a sequence of computations with 4/h? fixed we have 

| Zit O(h?\ 


It is apparent that Eqs. (5.4) are of the same form as (3.1), with nonlinear and 
truncation error terms added on. Hence in the analysis of (5.4) we only need to 
consider the effect of these extra terms, since # and v in (3.1) are replaced in (5.4) 
by U and V, which are a fortiori Lipschitz continuous. 

We intend to prove by induction that the f,-norm of the error [\|w’"|? + 
||u’” |? + |v’”|?}# is O(A?). This is done by showing that (5.4) is stable in f, under 
the induction hypothesis. We must estimate inner products of the form 


(o£ o9!* 3, AAge (a! ™Q" + u'r") + 2Agy (v'™Q" + 0"%0'") + ALT,). 


The errors uw’ and v’ are estimated as follows: 


Defining 
Do,=hAg,, Doyy=hAg,y, 
from 
U+u'=D,(P+y); V+ =—Dy,(P+y’) 
we obtain 


u! =Dy,y’ +0 (h): v' = —Do, py’ +0(h). 


Let « be the minimum eigenvalue of the difference operator —V? in R. Its 
magnitude will depend on the size and shape of the obstacle, but will be quite 
large for a rectangular obstacle whose size and position are commensurable with 
the width of the channel. The Rayleigh quotient bound 
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leads to the bound for yp’: 


Psat’, —V?y') = —a4(y’, T,—o') Sa4|y'|| 


yy 1b —«'|. 


Now take the inner product of (5.4b) with y’”, apply (A.5) and we have 


[Div "P+ Diy" P= — (py, T, —0'") 
Sy" IZ, -o""| (5.5) 
Sa? (|Z, + lo’? 
Thus ||v’”|| and |v’"|| are O(h2) by hypothesis. 
By the inequality between the ¢,, and the /,-norms 
|e *|, |e’ */ =O). (5.6) 


Hence w and v are Lipschitz continuous and the assumption we made in § 5 is 
justified. 
The additional inner products are expanded using the rule 4),(u’Q)= 
Ay, 2+,’ +y,Q-Ao,u', etc. to obtain 
(@’ "474+ @'"?, Bue ted eel, OF uo.” 

+ (U,2" — 2") Aye! + (1,.2" — 2") Ay, 0! 

+92" (A9,4'" +A, 0"”) 

+ Ags (u’"e!*) + Ag, (v'*™) 

+AtT,). 
Since 
0U oV 
Ox st oy 


Ao, @" + Ay, v” =0= 
we can show by Taylor expansion that 
Aor ™ ” + Agyv'” =AyyU” + Ay, V"® =O (h*). 
Hence on applying (5.6) we bound the inner product by 
a fe + fo 
of + Jo" + 0H4)). 


YA 


One 


AO (h) []o’""*P + lo" P+ | 
+ [20 (h) +0(A2)] ( 


Define 
uw’ n 


E,,=|o’"P +o" P+ lo" P+ [oP 
+ @?D, + 24(1+ 9) (@’", Ao,(Uw')" + Ao, (Vo')") 
V\)< d, and E,, =0 (h') 


, 


where Dj, has the same form as D,,, in w’. If A max(|U 
we can show that E,, is equivalent to the /, norm and 


E,, =O (Ab) (lo’* |? + lo" P + |e" P + [lo P) 
+0 (At) (llo’"**/P + jo’? +0 (74)) 
=0(At) (Exy1 +E,) +AtO (4). 


E 


n+l 
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Hence constants L and M exist such that for sufficiently small A?, 


Sats Ee Te, EIA 
and thus 
E, Se? (Ey) +14 Mh] 
=O(h*) forallmsuchthat nAtST, 
given that the initial error Ey is O(h4); i.e. O (h?) in f,. Our induction assumption 


is thus justified, and the errors about a C* solution of (2.1) are O(h?) in f,, and 
hence O(h) pointwise in the interior of R. 


Appendix 
Let uw and v be functions defined on a mesh with spacing h in [0, L] x [0, 1] 
such that u,;=«(th, 7h) for OSt SIT =L/h, 0Sj7 SJ =1/h. The following rela- 
tions can then be proved: 
In the x-direction we have 


(4, Mor?) = (Mor, V) + h/2 (Mya 5015 — M75 U1-1,; — 405%; +%1j%;)  (A-1) 


with the corresponding result in the y-direction. 


Next, 
(4, Ag?) = —(Ag,%, ¥) + hf2 (upg, , 075 + M15 Vz-1, 7 — 405215 + %1j5%;)- (A-2) 
When w and v are periodic in x, (wu, 4),v) = —(Ay,u, v) and in particular 


(4, Ap, %) =0. (A.3) 
Also in the periodic case ||T;,~|?= |||? and hence 
Ao. PS2(\Z, «|? +|Z_, uf) = |uP. (A.4) 
If uw is constant at the boundaries 
(u, 0,4) = —||4,, uP. (A.5) 
In two dimensions we have 
(Mx T byt, Agx + Agy %) =a. (A.6) 


These identities are modified when the inner product involves a variable 
coefficient a, which we assume to be Lipschitz continuous. Using the identities 


Ao, (au) =aAg,u+#4(T,,u4,,a+T_,uA_,a) 
Hox (4M) = 4 po, t +3(T,,wA,,a—T_,wA_,a) 
we find that (A.2) generalises to 
(4, Ags (@2)) = — (% Aox(au)) +0 () |} >| 
+ h[2( toy, 5 (40) 75 + ty ; (20) 7-4, ; — U0; (42), ; — Hj (42) 0) 
and (A.6) becomes 


(A.7) 


(u, + My 4, Ag, + Ap, (au)) =O (h) ule. (A.8) 
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Nichtnegativitats- 
und substanzerhaltende Differenzenschemata 
fiir lineare Diffusionsgleichungen* 


RUDOLF GORENFLO 


Eingegangen am 25. Juli 1969 


Abstract. Explicit difference schemes preserving non-negativity and mass as 
does the described diffusion process are given for the Fokker-Planck partial differential 
equation (Kolmogorov’s forward differential equation) in IR”. These schemes can be 
interpreted as descriptions of discrete non-stationary inhomogeneous random walks. 
They are stable in the maximum norm, and if an initial condition is given, they 
allow approximate solution of the differential equation as well as approximate Monte- 
Carlo simulation of the basic diffusion process. Essential conditions are the bounded- 
ness of the coefficients of the differential equation and of some of their derivatives 
and a strong diagonal dominance of the diffusion matrix. The results can be generalized 
to general linear parabolic differential equations (conservation of mass being lost, 
of course). By means of the developed method it is also possible to construct difference 
approximations of non-negative type to linear elliptic operators in IR” if these operators 
satisfy analogous conditions. 


1. Einleitung 


Diese Arbeit behandelt nichtnegativitats- und substanzerhaltende explizite 
Differenzenschemata fiir die allgemeine Fokker-Planck-Gleichung im IR”, »=1, 


De Ee @ a 
i Say 2 Bx; Ox, (Pin™) =o ax, 41%) 
j= = j= 


und fiir die Fokker-Planck-Gleichung mit selbstadjungiertem elliptischem Anteil? 


(FPS) ieee Dies 


In beiden Fallen sei eine Anfangsbedingung 


(A) u(x Oe es) CIR 


gegeben, und man interessiert sich fiir die Lésung u(x, ¢) bzw. eine Naherungs- 
losung in D =R" x {t|O St ST} mit einem festen T>0. Ein Punkt des R” mit 
den Koordinaten x,, %,,..., x, wird der Kiirze halber mit x bezeichnet. Im 

* Diese Arbeit entstand im Rahmen des Vertrags zwischen dem Institut fiir 


Plasmaphysik (GmbH) und der Europadischen Atomgemeinschaft iiber die Zusammen- 
arbeit auf dem bie der Plasmaphysik. 


ab 
1 Mit a;= > se 7 geht (FPS) in (FP) iiber. 
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folgenden sei stets (x, t)€D. Alle auftretenden Variablen und Funktionen seien 
reell. 

Die Matrix B(x, t) = {b;,(x, t)} set symmetrisch (also bj, =, ;) und gleichmafig 
positiv definit, so daB (FP) und (FPS) gleichmiéBig parabolische Differentialglei- 
chungen sind. B ist gleichmaBig positiv definit, wenn der kleinste Eigenwert von 
B nie kleiner wird als eine feste positive Zahl und wenn alle ,, beschrankt sind. 


Wir miissen noch das Erfiilltsein gewisser Glattheitsbedingungen fiir die Koef- 
fizienten 6;,(x, t), a;(x, t) und die Funktion g(x) fordern, damit die Lésung w(x, #) 
eindeutig existiert und ebenfalls so glatt ist, wie es fiir die in den nachsten Kapiteln 
durchzufiihrenden Konsistenz- und Stabilitatsiiberlegungen notwendig ist. Aus- 
sagen tiber die eindeutige Existenz der L6sung und den Zusammenhang der Glatt- 
heit der Koeffizienten und der Anfangsbedingung mit der Glattheit der Lésung 
findet man in [6]. 

Fiir unsere Abschdtzungen bendtigen wir die Existenz und Beschrankthett der 
b;, und gewisser threr Ableitungen nach x bis zur Ordnung 4, der a;, und gewisser 
threr Ableitungen nach x bis zur Ordnung 3, der Losung u und gewisser threr Ab- 
lettungen nach x bis zur Ordnung 4 und threr Ableitungen nach t bis zur Ordnung 2 
(hierzu miissen natiirlich auch die b;, und die a, beztighch t geniigend glatt sein). 

Wir verzichten darauf, moglichst schwache V oraussetzungen zu formulieren, damit 
diese Erfordernisse auch fiir die Lésung wu erfiillt sind, sondern begniigen uns mit 
der Feststellung, daB sie erfiillbar sind. (Man vgl. Kapitel 3 des Buches [6].) 
Wir nehmen sie im folgenden stets als erfiillt an, ohne dies immer ausdriicklich zu 
erwahnen. Wir verwenden gelegentlich den etwas vagen Ausdruck ,,hinreichend glatt, 
der die Existenz, Beschranktheit und Stetigkeit der jewetls bendtigten Ablertungen 
postulieren soll. 

Speztell g(x) set beschrankt und habe beschrankte Ableitungen nach x bis zur 
Ordnung 4. Wenn man dann fordert, daB u(x, t) beschrankt ist, so sind die Cauchy- 
schen Anfangswertprobleme (FP, A) und (FPS, A) in der Maximum-Norm korrekt 
gestellt. 

(FP) und (FPS) beschreiben die Diffusion einer im JR” verteilten Substanz, 
wenn man (x, t) als die Dichte dieser Substanz am Orte x zur Zeit ¢ auffaBt. Es 
bietet sich deshalb der Versuch an, zwei wesentliche Eigenschaften dieses Dif- 
fusionsprozesses, namlich die Nichtnegativitatserhaltung? und die Substanz- 
erhaltung? durch ein explizites Differenzenschema zu imitieren. Dies fiihrt auf 
Differenzenapproximationen nichtnegativen Typs, die auch unabhangig von der 
angedeuteten physikalischen Analogie von Interesse sind. Man erhalt namlich 
gleichzeitig Approximationen nichtnegativen Typs fiir die auf der rechten Seite 
von (FP) und (FPS) stehenden elliptischen Differentialausdrticke. Solche Ap- 
proximationen sind niitzlich bei der numerischen Behandlung nichtlinearer ellipti- 
scher Differentialgleichungen und bei Monte-Carlo-Methoden fiir lineare elliptische 
Randwertprobleme (man vgl. [1, 2, 4, 8, 12]). 


2 Nichtnegativitatserhaltung bedeutet: Ist u(x, s) 20, so ist auch u(x, ¢) 20 fiir 
t>s. Man vgl. [6], S. 42ff. Substanzerhaltung bedeutet: Die Substanzanderung in 
einem Gebiet V¢ R” mit glattem Rand kommt dadurch zustande, da® tiber diesen 
Rand Substanz zuflieBt oder abflieBt. Die Substanzerhaltungseigenschaft lat sich 
leicht mittels des Gau8schen Integralsatzes einsehen. 
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FaBt man w(x, t) als Wahrscheinlichkeitsdichte des Aufenthalts eines irrenden 
Teilchens am Orte x zur Zeit t auf, so beschreiben (FP) und (FPS) stetige Markov- 
Prozesse*® (man vgl. [3]). Die gegebenenfalls existierenden nichtnegativitats- und 
substanzerhaltenden Differenzenapproximationen kénnen dann interpretiert wer- 
den als Beschreibung diskreter Irrfahrten, und man kann sie Differenzenschemata 
vom Irrfahrttypus nennen ([7]). Sie gestatten also neben der approximativen 
numerischen Lésung des Anfangswertproblems auch die approximative Monte- 
Carlo-Simulation des Diffusionsprozesses. Im konkreten Einzelfall kénnen natir- 
lich beziiglich Ausnutzung eines Computers effektivere Methoden vorhanden sein *. 
Im Falle n =1 sind einige einfache Beispiele diskreter Irrfahrten und ihr Ubergang 
bei verschwindender Schrittweite in Diffusionsgleichungen schon lange bekannt 
(man vel. z.B. [5], S.354ff., [41, 13]). Fir den Fall »=2 mit nur von ¢ ab- 
hangenden 6,, und a; vgl. man [9]. 

In diesem Zusammenhang sei auf die wahrscheinlichkeitstheoretische Bedeu- 
tung von B und a hingewiesen. Bezeichnet x = x(¢) die Bahn eines gema8 dem 
durch (FP) beschriebenen Markov-ProzeB irrenden Teilchens, und ist <-> der 
Operator der mathematischen Erwartung, ist ferner Af>0 und Ax=~x (t+At) —x(?), 
so ist fiir At-+0 

<Ax> = a(x, t) At+o0(A?t), 
<Ax;Ax,> = b;,(%, t) At +0(A?2). 
wy 
Man nennt a= | : | den Driftvektor, B= (b,,) die Diffusionsmatrix. 


ay, 


Theorem 2 von Motzkin und Wasow ({12]) la8t hoffen, unter ziemlich allge- 
meinen Bedingungen Differenzenschemata der gewiinschten Art finden zu koénnen; 
beschrankt man sich aber auf ein von vornherein gegebenes System von Nachbar- 
gitterpunkten, so lehrt Theorem 1 von [12], daB man ziemlich einschrankende 
Bedingungen an die elliptischen Differentialausdriicke auf den rechten Seiten 
stellen muB. In dieser Arbeit sollen Schemata einer noch zu prazisierenden ein- 
fachen Struktur — und hinreichende Bedingungen fiir ihre Existenz — angegeben 
werden. Wir werden uns beschranken auf kubische Gitternetze im IR” mit den 
Gitterpunkten x=mh, h>0, m ein Multiindex, und auf das Nachbarschafts- 
system von %, das in diesem Netz alle Punkte mit den Abstanden 0, # und \2h 


3 Man nennt dann (FP) auch die Vorwartsdifferentialgleichung von Kolmogorov. 

4 Bekanntlich ist ja auch die Nichtnegativitat der Koeffizienten eines expliziten 
Schemas fiir parabolische Differentialgleichungen keine notwendige Bedingung fiir 
Stabilitat ([4], S. 113ff.). Man wird auch, falls nur ein endliches Gebiet vorliegt und 
Randbedingungen gegeben sind, mit Vorteil implizite Methoden benutzen. — Eigent- 
licher Anla8 fiir die Entstehung dieser Arbeit war eine Vortragsreihe iiber Diffusions- 
prozesse, die der Verfasser vor Physikern des Instituts fiir Plasmaphysik hielt. Er 
beschaftigte sich in diesem Zusammenhang mit approximierenden diskreten Irrfahrt- 
prozessen und stellte sich die Frage, wie weit man die in [5] fiir »=1 beschriebenen 
Irrfahrt-Approximationen verallgemeinern kann. Da er sich gleichzeitig in die Theorie 
der Differenzenschemata einarbeitete, benutzte er den Apparat dieser Theorie zur 
Untersuchung von Konvergenz und Stabilitét. Er hofft, den Nutzen der beschriebenen 
Schemata bei numerischen Rechnungen und zweckmaBige Modifikationen fiir para- 
bolische Randwertprobleme in einer spateren Arbeit untersuchen und diskutieren zu 
k6nnen. 
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von x enthalt. Als wesentliche Bedingung fiir die Existenz entsprechender Dif- 
ferenzenschemata wird sich die Diagonaldominanz der Diffusionsmatrix B er- 
geben. Mit dieser Bedingung ist eine wichtige Klasse von in den Anwendungen 
wichtigen linearen Diffusionsgleichungen erfaBt. 


2. Konsistenz und Stabilitat des Differenzenschemas, 
auch fiir allgemeinere lineare parabolische Differentialgleichungen 
Wir geben der Vollstandigkeit halber einen kurzen Abri8 (in Anlehnung an 
[10], Kapitel 9) der Theorie von Konsistenz, Stabilitat und Konvergenz fiir das 
inhomogene Anfangswertproblem 


(1) MS Aee Y ese Sy, 


“u=4«4(x,#) gesuchtin 0OS/ST. 
Dieser AbriB ist zugeschnitten auf unsere Bediirfnisse, enthalt aber im wesent- 
lichen wohlbekannte Aussagen. A= A(x, ¢) sei ein linearer Differentialoperator, 
der nur Ableitungen nach den Komponenten *; von x enthalte (auch Ableitungen 
hoherer Ordnung sind zulassig). Das Problem (1) sei korrekt gestellt in der 
Maximum-Norm, d.h., w soll eindeutig existieren und normstetig von den Anfangs- 
werten abhangen — in einem geeigneten Funktionenraum. Fiir irgendeine be- 
schrankte Funktion w(x, ¢) ist die Norm so definiert: 


|o (¢)|| = sup |w(x, 2)|. 
x€ IR" 


Wir diskretisieren so: Zu h>0O und t>0 mit tX)(h), AShy, T(h) eine 
geeignete Funktion von /, betrachten wir die Punkte x =mh, t=y t des Gitters 
D =D (h, t). Dabei ist m ein Multiindex mit den » ganzzahligen Komponenten m,, 
—co<m;<00,7=1,2,...,%, und y durchlauft die ganzen Zahlen aus OS yST/t. 
Es interessiert der Grenziibergang h— +0. 


Setzt man L = ae —A, so geht (4) iiber in 
(2) Lu=/f, u(x, 0)=¢(x))  « gesuchtin 0OStsT. 
Konsistenz einer Differenzenapproximation mit dem Gitter D 
(2’) LypU=F, U(x,0)=G(x), U gesucht in D=D(h, t) 


zu (2) bedeutet nun, daB fiir die lokalen Abbruchfehler® 4, , 7 in 


(3) Lpu=f{+a, F=f+o, G=gt?7 

bei h—->0 gilt 

(4) sup |A([p>>0, sup |@o>0, sup |7Olb>0. 
OS/tST 0stsT 0stST 


5 Es ist zweckmaBig, hier Funktionen F und G einzufiihren, da es im Sinne der 
diskreten Imitation der zugrunde liegenden Diffusionsprozesse liegt, anstelle von / 
und g Mittelwerte iiber gewisse mit dem Gitter zusammenhangende Wiirfelbereiche 


zu nehmen. 
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Man kann auch so sagen: Der Differenzenoperator Lp ist eine O (h*)-konsistente 
Approximation an den Differentialoperator L, wenn fiir hinreichend glatte Funk- 
tionen v(x, ¢) gilt 

(Zp —L) up =O). 


Dabei sei p> 0. Gewiinscht ist, daB bei A-+0 die Lésung U von (2’) gegen die 
Lésung wu von (2) strebt. Wir miissen noch die Gitternorm |- ||, definieren. Fiir 
irgendeine beschrankte Funktion w (x, ¢) ist 
l2@lo=le@lba,.= sup |w(x, 4). 
xER"™ DD (h, Tt) 

Nun wird in der wirklichen Rechenpraxis nicht U aus (2’) ausgerechnet, 
sondern es wird, mit lokalen Rundungsfehlern g und y, eine Gitterfunktion W 
aus dem Problem 


(27) LpW=F4+q, W(x,0)=G(x)+y,  W gesucht in D, 
ausgerechnet. Es gilt das 


Stabilitatslemma 1. Wenn fiir beliebige beschrankte Funktionen F und G aus 
(2') folgt, daB mit positiven Konstanten K,, K,, die von den Inhomogemtaten F, G 
und von der Feinheit des Gitters D= Dh, t) unabhangig sind, 


(5) IPOb SK, sup [FOI +KalGlo,  ¢=r7, 


gilt, so gilt fiir die Losung W des Problems (2") die Ungleichung 
(6) I4%@) = 4Olb= 5, sup [b-> — Aba Kale typ. 


Das Lemma lehrt, daB bei rundungsfehlerfreier Rechnung W(x, ¢) gleichmaBig 
gegen u(x,t) konvergieren wiirde, und da8 die auftretenden Rundungsfehler, 
deren GréBenordnung man in der Praxis oft schatzen kann, nicht beliebig stark 
vergroBert werden, wenn (5) gilt. (5) ist die definierende Bedingung fiir Stabilitat 
des Verfahrens. 

Zum Beweis des Lemmas braucht man nur zu beachten, daB Lyj(W —u)= 
g+g—A und W(x, 0) —u(x, 0) =p +y ist. 

Wir formulieren noch als hinreichende Bedingung fiir die Giiltigkeit von (5) 
das bekannte ,,Indexkriterium“ als 


Stabilitatslemma 2. Wenn mit einer von der Feinheit des Gitters unabhingigen 
Konstanten C =0 aus (2') 


(7) |U@+)bSU+Cr) UO) +7|FO| 
folgt, so gilt (5) mat 


ID» t— VE, 


K,=(e7—1)/C fir C>0, K,=T fir C=0, . K,=€". 


Mit etwas handlicherer Bezeichnungsweise folgt dieses Lemma sofort aus fol- 
gender Aussage: Es seien {s,} und {v,} Folgen nichtnegativer Zahlen mit v, < 


(1+-C t) up 8,24: Ferner sei. C 20,.N = [7/2] 0S95 NPS max siDann 
O0SrvSN-1 


ist v, SKyvy + K,S. Diese letzte Ungleichung ergibt sich so: Eine einfache Re- 
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kursion und anschlieBende Abschatzung liefert 


y—1 
v= (44 Cr)", + 7S Dd (44 Cr) 


j=1 


N . 
<(ic-Or) nS (tC cr) (=.Ko0,4- Ah. - 
j=1 
Im folgenden werden wir das hinreichende Stabilitatslemma 2 benutzen. Jetzt 
wenden wir uns unseren Diffusionsgleichungen zu. (FP) und (FPS) lassen sich 
durch Ausdifferenzieren umformen in eine lineare parabolische Differentialglei- 
chung der Gestalt 


me 
ou 


4 ” Ou “, OU be 
(8) “a EM mit Pu=— pap? Des 5 +éu. 


Fiir die Diskretisierung D sei, mit spater noch weiter einzuschrankendem kon- 
stantem o> 0, 
(9) T=26n". 


Das Differenzenschema, nach U(x, t+ 7) aufgeldst, lautet 


Ue t= ty > oe tk) Orr h.t); t=, —v= 11h, 
(10) reR 
y=0,1,2,..., [T/t] —1. 


Dabei sei i ein festes System von Multiindices 7 mit 7’= (1, 72, ..., 7). Setzt 
man, um den AnschluB an die weiter oben dargestellte Konsistenz- und Stabilitats- 


0 


theorie zu finden, d= P und L= ~ —P, so ist fiir eine hinreichend glatte 
Funktion v(x, ¢) sinngemaB 


1 Rue 
Lpv=— {u(x,t+7) —v(x, oe Pr thier) \o(x+rh,t) +28 


Wegen 
= ~ {o( %,¢ +t) —u(x,t)} ==, +0(t) 
ergibt sich als 0 (41)-Konsistenzbedingung 
(Pp — P)v=o0(1). 


Da diese fiir jede hinreichend glatte Funktion v erfiillt sein muB, ergibt sich 
durch Taylor-Entwicklung ihre Aquivalenz mit folgenden 3 Bedingungen (dabei 


stip, —= p(s, th 0): 


(K 1) e(de—1)>e 
rER 
ji 
(K 2) = De D: >a; >) bei h>0, 
reR 
he 
(K 3) ayalaei de —> 5, 
A 


318 


454 R. Gorenflo: 


bzw. 
),=1+Et+o(t), Dr; p,=20ha;+o(h), 
reR reR 
me 1; V_ Py = 26 b;, +0(1). 
reR 


Notwendig und hinreichend fiir Nichtnegativitats-Erhaltung ist das Erfiilltsein 
der Bedingung 
(NN) BAt, ti Nh, tT) a0, Teer, 


die besagt, daB die Approximation von nichtnegativem Typ ist. 

Wir haben jetzt die Mittel bereitgestellt, einen allgemeinen fast trivialen Satz 
iiber lineare parabolische Differentialgleichungen und ihre Approximation durch 
Differenzenschemata nichtnegativen Typs auszusprechen. 


Satz 1. Das lineare gleichmafBig parabolische Anfangswertproblem 


(8’) ot = Putf mit P wiein (8), u(x, 0)=—g(x), 


hinreichend glatten Koeffizienten, hinreichend glatter Funktion f =f (x,t) und hin- 
reichend glatter Anfangsbedingung g(x) werde approximiert durch ein Schema 
(10’) U(x,t+1) => ?,(x, t; h, t) U(x +rh, t) +4F (x, t) 
reR 

mit U(x, 0) =G(x). Gesucht ist eine Naherungslosung in 0 St ST. F und G seien 
konsistente Approximationen an f und g, (10’) erfiille die Konsistenzbedingungen 
(K1), (K2), (K3) und die Nichtnegativitats-Bedingung (NN). Dann ist (10’) ein 
stabiles Schema®, das zu (8') konsistent tst. 

Beweis. Aus (10’), (8’) und (NN) folgt 


|U¢+[b>S(¥ A) [OOo +7 IF Olo- 


Mit (K1) ergibt sich nun sofort (7). Zu ¢>0 gibt es ein hg >0 so, daB man 
C = sup |¢| +e fiir 4<h, nehmen kann. 

O<tST 

In Satz 1 ist es nicht erforderlich, daB die verkiirzte Gleichung du/ét= Pu 
eine Fokker-Planck-Gleichung ist (B muB8 natiirlich symmetrisch und gleich- 
maBig positiv definit sein). Die allgemeine Gl. (8’) entsteht, wie leicht nachzu- 
rechnen ist, aus einer Gl. (FP) durch additive Hinzufiigung eines Terms c (x, t) u 
und eines Quellterms /(x,¢). Die Interpretation dieser Terme als Entstehungs- 
dichte von Substanz liegt auf der Hand. cu bedeutet Entstehungsdichte pro- 
portional zur gerade vorhandenen Substanzdichte. (Ein negativer Term c u oder 
7 entspricht einer Vernichtungsdichte von Substanz von der GréBe seines Be- 
trags.) Hat man zu (FP) eine Approximation, in der 9(x, t; h, t) nicht kleiner 
als eine positive Konstante ist (dies werden wir spater durch hinreichend kleine 
Wahl von o erzwingen kénnen), so ist fiir hinreichend kleines # (und damit hin- 
reichend kleines t) auch f)+-c t > 0, die Approximation bleibt also von positivem 
Typ, wenn sie es vorher schon war. Aus diesem Grunde geniigt es, sich auf 


6 Beziiglich des Falles einer Raumdimension (7=1) vgl. man [4], S. 107 ff. 
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Gln. (FP) zu beschranken (von denen (FPS) einen Sonderfall darstellt). Man 
erhalt durch auf der Hand liegende Modifikationen der darzustellenden Theorie 
auch Aussagen tiber Approximationen nichtnegativen Typs fiir allgemeine in- 
homogene lineare parabolische Differentialgleichungen (8’). 

Wir prazisieren den im Gitter D =D (h, t) zu verwendenden Differenzenstern, 
d.h., das zu einem Punkt « im IR” gehdrende System von Nachbargitterpunkten? 
N(x) = N(x; h). Der Einfachheit halber sei {9(«)} symmetrisch und translations- 
invariant, also 
(11a) mit xEeX(y) sei yER(x), 


(11b) M(x) = MO) + x. 
Der Zusammenhang mit dem Indexsystem ft ist gegeben durch 
R= {r| x +rheN(x, h)}. 
Statt (10) schreiben wir, mit der Bezeichnungsweise’ 
2=—£-+-9h, p(x, 8) ht) = f(s, X; 1) — P28 th, zc) 


U(x, t4-2)—= Dd ple, a7 Ul2e,0); tre, 
(12) ZEN (x) 
Per O 4, 2, flat hilt 


Die Nichtnegativitatsbedingung (NN) geht iiber in 
(13) o(%,y;t)20 fir yedt(x). 


Aus formalen Griinden ist es zweckmabig, p(x, y; 4) =0 zu setzen fiir y¢N(x). 
Man kann dann in (12) und (14) unterm Summenzeichen ¥t(x) durch R’,7D 
ersetzen. 

Denkt man sich fh” U(z,¢t) als eine zur Zeit ¢ im Gitterpunkt z befindliche 
Substanzmenge und #(z, x; ¢) als den Bruchteil dieser Substanzmenge, der im 
Zeitintervall (¢,¢-+-7) von z nach x wandert, so ist Substanzerhaltung gleich- 
bedeutend mit der Bedingung 
(14) D P(*, y3 ) =1. 


yEN(x) 
Im Sinne dieser Auffassung liegt es, in der Anfangsbedingung 
(15) U(x, 0) =G(x), 


x+oh 
(16) Gojeay fread te mit (ov (54, .7) 1/2 
x—oh 
zu nehmen. Dies ist O (h?)-konsistent mit u(x, 0) =g(x), da 
n k—-1 


n 
6) el s(&ER, +S 3 Ra) 
j=1 k=1 j=1 
7 Aus Griinden einfacherer Schreibweise sei in Zukunft die Abhangigkeit des 
Nachbarschaftssystems und der Koeffizienten p(z, x; t) von h und t normalerweise 
nicht angeschrieben. 
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ist mit a2g(x) 
OX; OX, 


R 


jk— 


raver er 
Die Aufgabe ist nun die, (13) und (14) erfiillende konsistente Schemata zu 
finden. Hierzu beschriinken wir uns von jetzt ab aut einen einfachen Typ des Gutters: 
MN (0) =M(0; A) bestehe aus den 1+2n® Punkten 0, the;, Ahe;+phe, mit 
j+hk, wobei (A, ) die Zahlenpaare (—1, —1), (—1,1), (4, —1), (1,1) durchlauft. 
e, ist der Einheitsvektor in x,-Richtung; zwecks Abktirzung der Schrevbarbeit ser im 
folgenden h;=he,. 


3. Diskretisierung der Fokker-Planck-Gleichung 


Die unmittelbar sich anbietende, O (h?)-konsistente Approximation von (FP) 
hat die Gestalt 


— {U(x, t+) — U(x, 9} 


Ki 2, {b;;(x +h,;, t) U(x +h,, t) 
j= 


(17) — 2b; ;(%, t) U(x, t) + b;;(%¥ —h,, t) U(x —h,, t)} 
nm j—1 
co te Dah Die dale + Ah; + phy, t) U(% + Ah; + whp, b) 
J=1 k=1 Ap 


Sh Ds (ai (* +h; 0) U(x +h,, t) —a;(x —h;, t) U(x —h;, Oy. 


Hier und im folgenden bedeute >’ stets Summation iiber die am Ende von Kapitel 2 
4, 
angegebenen 4 Zahlenpaare (A, mu). : 

Der Diskretisierungsfehler ist gleich h? multipliziert mit einer gewichteten 
Summe aus vierten (auch gemischten) rdéumlichen Ableitungen der 6;,u, dritten 
radumlichen Ableitungen der a;u an Zwischenstellen (&, ¢), und der zweiten zeit- 
lichen Ableitung von u an einer Zwischenstelle (x, s). Die Stellen & liegen in der 
konvexen Hiille von (x). Der Vollstandigkeit halber sei hier erinnert an 


(9) t= 20h'. 
Aus (9) und (17) ergibt sich 
(18) Ula, 1) =>) az, a: f)U(z.2) 
ZEN (x) 
mit 
(19a) q(x, x; t)=1—20 D).b,, (x, 0), 
j=1 
(19b) q(x, x h;; t) =o{b;;(*, t) tha,(x, dh, 
(19c) q(x, x + Ah; + why; t) = Apdj, (x,t), 7 ek 


Aus (17) findet man die GréBen ¢(z, x; t). Wegen der Symmetrie-Eigenschaft (11 a) 
des Nachbarschaftssystems bereitet es aber keine Schwierigkeit, sich zu iiber- 
legen, wie die GroBen g(x, z; t) in (19) aussehen. 
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Die q(x, y; ¢) erfiillen die Substanzerhaltungsbedingung; man rechnet leicht 
nach, daB 


(14’) 2X 9(% y3 t)=1 

VEN (x) 
ist. Die Nichtnegativitatsbedingung jedoch ist i.allg. leider nicht erfiillt, da unter 
den GroBen (19c) negative auftreten, wenn nicht alle 0,, fiir 7+ verschwinden 
(Aw nimmt ja die Werte +1 und —1 an). Um diesem Ubelstand abzuhelfen, 


suchen wir Korrekturen c(x, y; f) und Bedingungen ihrer Existenz derart, daB 
die GréBen 


(20) P(x, vs = a(% y5 t) Fe(x, y; 2) 
den Bedingungen (13) und (14) geniigen. Es soll also gelten 


(13’) q(x, y; t) +c(x, y; t) 20, 
(14’) pe(envs tO. 
y EN (x) 


Ferner wiinschen wir, daB dann (12) ebenso wie (18) zu (FP) O (h?)-konsistent ist. 
Subtrahiert man auf beiden Seiten von (12) U(x, ¢) und dividiert man dann durch 
T, so ergibt sich wegen (9) und (20), daB fiir O (h?)-Konsistenz die Bedingung 


(21) » c(z, x; t) v(z, t)=O(A4), v(x, ¢) hinreichend glatt, 
zEMN(x) 


notwendig und hinreichend ist. 


Die Grundidee fiir die Anbringung der Korrekturen ist die, zuerst anstelle 
von (19c) nichtnegative GréBen einzufiihren. In (19b) und (19a) miissen dann 
die in (19c) angebrachten Korrekturen wieder kompensiert werden, und zwar 
gerade so, daB (21) gilt. Nichtnegativitat der korrigierten GroBen (19b) und (19a) 
erreicht man dann bei hinreichend starkem Uberwiegen der Diagonalelemente 
b;;(x, ¢) tiber die Betrage der anderen Matrixelemente durch Wahl hinreichend 
kleiner / und o. 

Wir wahlen fiir 7k Funktionen g,,(x,¢) mit der Symmetrieeigenschaft 
gia (%, t) =, ;(x, t), die zusammen mit ihren als existent vorausgesetzten Ablei- 
tungen nach x bis zur Ordnung 4 beschrankt seien. Dann erfiillen die Korrekturen 


(22a) C4: ot) 03 (5+ D ) gale t), 
Gis 

(22b) o(x, «th —o(5 120y ) ils ‘), 
k=j+1 

(22c) e(x, x +h, +uh,; ) = ets ), 7k, 


die Bedingung (14’) und, wie eine langere elementare Rechnung zeigt, auch die 
Bedingung (21). Die in (21) links stehende Summe ist gleich 44 multipliziert mit 
einem Ausdruck, der Ableitungen der Funktionen g,,(x,¢) und u(x, ?¢) nach x 
bis zur Ordnung 4 enthalt, genommen an geeigneten Zwischenstellen der kon- 
vexen Hiille von 9t(x). 


458 R. Gorenflo: 


Damit fiir 7 +k die p(x, x +Ah,;+h,; t) 20 sind, fordern wir das Bestehen 
von 
(23) gin (%, t) ZO, (%O|, JR. 


Um p(x, *+h;;t)20 zu haben, muB, da i.allg. nicht alle a, verschwinden, 
jedenfalls b;;> 2 Sik sein. Man kommt zum Ziel, wenn man mit Konstanten 


6>0,¥, ipaiea 1, nd einer Funktion 7 (x, t) =6 das Erfiilltsein der Bedingungen 


t : 
(24) [Bia (2, 8) Saja (ts) S[Bn(x, | +22, j+k, 


(25) b;;(x, 1) = ba +3) |;,(%,2)|+4(*,t), 7=1,2,...,%, 


fordert. Die Diagonaldominanz (25) ist die entscheidende einschrankende Be- 
dingung fiir die Diffusionsmatrix B(x, ¢). Aus ihr folgt nach dem Satz von 
Gerschgorin ([14], Kapitel 1), daB der kleinste Eigenwert von B (x, t) stets > 0 ist. 


Aus (24) und (25) folgt 
P(x, x+h,; t) 2o{(1 —y) 6 —hla;]} 


und dies ist =O, wenn 


(26) 0<hsSh=(1—y) é/ sup |a;(x, 2)| 


j, £,0<tST 


ist. Um p(x, x; ¢) =O zu erhalten, nehme man 


er) o<esm sup 5 f5,,(2.)-2(5 + > Jelena. 


*,0OStST j=1 k=j+1 


Der Ausdruck hinter sup >) ist nach (24) und (25) 


#,0StST j=1 
1» 1 (5 soidiny n 
Scacwry aabcullapeerans iow t pig WAM 
keh 
Mit Beriicksichtigung von Satz 1 ergibt sich als Resultat 


Satz 2. Es existiere eine Funktion (x,t) 2d>0 so, daB fiir 7 =1, 2,..., n die 
Diagonaldominanzbedingung (25) gilt. Man nehme fiir 7+:k hinreichend glatte® 
Funktionen g;,(x,t), &;,=8&,;, We den Bedingungen (24) geniigen. Fiir h und t 
sollen (26) und (27) erfiillt sein. Mit den Koeffizienten p (x, y;t)=q(x, y3t)+c(x, y;2), 
q aus (19a, b,c), ¢ aus (22a,b,c), ast dann {(12), (15), (16)} ein zw {(FP), (A)} 
O (h*)-konsistentes und stabiles Differenzenschema, das nichtnegativitats- und sub- 
stanzerhaltend wirkt. 


Anmerkung 1. Eine Stabilitatskonstante C des Schemas (man vel. das Stabili- 
tatslemma 2 in Kapitel 2) findet man durch Betragssummenabschatzung aus der 


Stabilitatssumme >) #(z, x; t) (man vgl. Satz 1). Es ist ja 
ZEN (x) 


(28) dD: P(2, 3) S14Cr. 


zEMN(x) 


8 Die g;, und ihre Ableitungen nach x bis zur Ordnung 4 sollen existieren und 
beschrankt sein. 
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Eine langere Rechnung gibt 


n n 
1 O* bj 0a 
b(z,*; 4) =44+75— = e981 j 
(29) ois 2 2 2: Oi; hy 2 Ox; 
n -—1 a £ 
Seas ( ee - Sp ae al STE 4 
4 k=1 Ox} OnR Ong 6 dx} 12 ox : 


Anmerkung 2. Fiir den lokalen Abbruchfehler des Verfahrens (man vgl. Ka- 
pitel 2) ergibt sich nach langerer Rechnung 


tT OK 


(30) (Lp —L) w=} “oa + (B+4,+Q). 


Mit Vin= jr, %=a,;u, B;,=b;,u ist hier 


n 2 n j—-1 
=. 1 O* Bj; 1 (=. 0 \4~ é 0 \4 x 
By 24 £1 ox} ia rp >{ Ox; 7 se Bin—(ae > aaa) Bis}. 


=1 k=1 


j=1 J 
qo j-1 OY; % OAD; 1 ta) 0 \4e 1 a a\4. 
: 7 me at ox} a, Ox} ai 2 G 25 aa) Vik et 2 & = Pal Fin} s 


Die Zeichen ~ bedeuten, daB die Ableitungen nach x an geeigneten Zwischen- 
stellen (nicht jedesmal an denselben Stellen) der konvexen Hiille von ¥t(x) zu 
nehmen sind, die Ableitung nach ¢ an einer Stelle (x, s) mit t<s<¢+t. 

Je nach Art der Anwendung des Mittelwertsatzes der Differentialrechnung 
kann man fiir Stabilitatssumme und Abbruchfehler natiirlich auch andere Aus- 
driicke erhalten. Wesentlich sind die Ordnungen der auftretenden Ableitungen, 
aus denen man erkennen kann, was unter ,,hinreichend glatt“‘ zu verstehen ist 
(man vgl. hierzu die in Kapitel 1 getroffene Vereinbarung). 


Sonderfalle. Wenn fiir 7=- alle b;,=0 sind, kann man g;,=0 und ii 
jnt (6 (0; ;/|@;|) nehmen. Wenn fiir ein Paar (j,k) mit 7--k entweder stets b;,=0 
tee stets b;,<0 ist, kann man das zugehirige 8; n= |5j;_| nehmen (effektiv ver- 
Rleinert sich dadurch der Differenzenstern um zwei Punkte). Wechselt keines der 
b;, sein Vorzeichen und sind alle a;=0, so kann man in (24) y=1, mm (25) n=0 
zulassen und auf (26) verzichten (aus (25) folgt dann allerdings nur noch, dag B 
positiv senndefinit rst). 

Es sei noch erwahnt, da8 bei Auffassung von (12) als Beschreibung eines 
diskreten Irrfahrtprozesses (h” U(x, t) = Wahrscheinlichkeit des Aufenthalts eines 
irrenden Teilchens am Orte x zur Zeit t, p(x, 2; t) = fasted) des Uber- 
gangs vom Orte x zum Orte z im Zeitintervall (¢,¢-+ 1), h” >)G(«)=1) mit 
At=vt exakt die Relationen #eRR" 


<Ax>=a(x,t)At, <Ax;Ax,> =);,(x, At 


gelten (man vgl. Kapitel 1). 
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4. Diskretisierung der Fokker-Planck-Gleichung 
mit selbstadjungiertem elliptischen Anteil 


Die Differentialgleichung (FPS) ist ein Sonderfall von (FP), sie geht mit 
ix Ob; k 
a, — 
1 2 reat Ox, 


in (FP) tiber. Man kann (FPS) also mit den Methoden des Ka- 


pitel 3 behandeln. Die spezielle Bauart von (FPS) empfiehlt aber ein anderes 
Vorgehen, bei dessen Darstellung wir uns aber kiirzer fassen konnen. 


Eine O (h?)-konsistente Approximation kann man so gewinnen: Man formt 
(FPS) um in 


ou 1 é ou 4wae a ou’ é é 
(31) a ~ @D > 0x; (0; aH | ap pes tee (a fal 0x; (Ge ok 


j=l 


und verwendet die Beziehungen 


oH =F fu(x,t4+7) —u(x, t)} +0(t), 


T 


5 (O43) = ge (Oise +S hi) (ule +h, 8) — aCe, 0) 


—bjy(x— shy, 7 (u(x, t) — u(x —h, i))} +0 (fh), 


é Ou 1 
= (Or a = ape (Dia (e + yj, t) (u(x +h + hy, f) —u(x thy —hy, f)) 
— bj, (% —hy, t) (u(x —h; +hy, t) —u(x —h; —hg, t))}+0 (2) 


fiir 7--k, und die Beziehung, die sich aus der letzten durch Vertauschung von 7 
und & ergibt. Ersetzt man wie tiblich w durch U, laBt die Terme O(t) und O (h?) 
weg und lést dann nach U(x, ¢+ 1) auf, so erhalt man eine Beziehung der Gestalt 
(18) mit 


(32a) q(x, %;t)=1 —0 2 {b;5(% + 5h;, 2) +0; ;(x —Eh,, D}, 
(32b) q(x, x+h,;; t) =0b,,(x+5h,,0, 
(32c) g(x, x +Ah; + phys t) = Apo (x + AK; t) HO (xtuhy, Dd}, G#R. 
Die Substanzerhaltungsbedingung >) g(x, y; 4) =1 ist erfiillt, im allgemeinen 
y EN (x) 


aber wieder nicht die Nichtnegativitatserhaltungsbedingung q(x, y; t) =0. Eine 
Besonderheit des Schemas ist die Symmetrie 


(33) q(x, ¥; t) =q(y, x; 2), 


die ein diskretes Analogon der Selbstadjungiertheit des elliptischen Anteils ist. 
Fiir die anzubringenden Korrekturen mu8 nun gelten 


(34) p(x, y; &) =q(x, y; t) +(x, y; 1) 20, 
(35) cn) =O. 


yEN(x) 
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ferner eine Konsistenzbedingung, die wieder die Gestalt (21) hat, wenn das re- 
sultierende Schema O(h?)-konsistent sein soll. Dariiber hinaus wiinschen wir 
Symmetrie 


(36) P(X. ¥; t) =O, X32) ow also. cle vs t)==c(y, x: 2). 


Man erreicht die gewiinschten Ziele, wenn man wieder an die Matrix B eine 
geeignete Diagonaldominanzbedingung stellt und wieder hinreichend glatte Funk- 
tionen g;, (x, ¢) mit g;, =g,, einfiihrt, die der Bedingung (23) geniigen. Man nehme 


(37a) (x, x: #) =) (= a > ) ial t)5 
j=1 R=1 k=j+1 
(37b) c(x,x+h,;t)= (5 +2 S |) teal (x,t) +ej,(e+h,, O}, 


(37c) AR Sco +Ah;+uh,; t) = = {ine Ay, t) + gyn (% +puh,, t}, 4=ER, 


Dann sind (36), (35) und (21) erfiillt. (34) ist erfiillt fiir die p(x, x +Ah,;+ph,; t), 
7 +k. Damit (34) fiir die p(x, x--/,; t) gilt, muB fiir alle x des Gitters und alle 7 
und fiirx—=-+1 


(38) > ee +3 | {gjn(% ) +e (4 +%h,, D)}S d(x + +— hy , 


sein. Aus Symmetriegriinden geniigt es, (38) fiir x= -+1 zu fordern. Addiert man 
die beiden Gln. (38) fiir x =-+1 und x=—1 und beachtet, daB g,;, 20 ist, so 
erhalt man 


(39) (S +S lente So, + bij.) £8j,(0— By 8. 
kajty 
Damit nun auch £(x, x; ¢) =O ist, braucht man nur noch o = 1/(2h?) hinreichend 
klein zu wahlen. 

(38) ist erfiillbar, wenn mit einer Funktion 7 (x, t) 20 die Diagonaldominanz- 
bedingung 
(40) ysl + By) S4(S +S \ (bale 0] + 1dpale +h OD mle O, 

i 2 


(ate ie PE 


erfiillt ist. Man nehme dann die g;, so, daB gilt 


7 ( 
(41) | bin (%, t)| = £5 \%, t) <|d,,(x, t) | “—- ae 
Als Resultat formulieren wir 


Satz 3. Es existiere eine Funktion (x,t) 20 so, dab die Diagonaldominanz- 
bedingung (40) gilt. Man nehme fiir 7 +:k hinreichend glatte Funktionen g;,(X, 1) , 
die den Bedingungen (41) und der Symmetriebedingung g;,—=&,; gentigen. Es sea 
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h>0 und 
Sel sup D Osle + By) +84, (0 — hj, t) 
»OStST j=1 
(42) 
ap 7 a ) gi (x, p}. 
b= bird 


Mit den Koeffizienten p(x, y; t) = q(x, y; t) +c(x, y; 1), g aus (32a—c), c aus 
(37a—c), ist dann {(12), (15), (16)} ein zu {(FPS), (A)} O(h?)-konsistentes und 
stabiles Differenzenschema, das nichtnegativitats- und substanzerhaltend wirkt. 


Anmerkung. Das Schema hat die Stabilitatskonstante C=O (man vgl. Ka- 
pitel 2). Wegen (35), (36) folgt nédmlich aus der Substanzerhaltung, daB die Sta- 
bilitatssumme >) £(z, x; t)==1 ist. Der lokale Abbruchfehler (Lp —L)wu ist h? 

ZEN (x, 
multipliziert mit ates Ausdruck, der von den Funktionen u, b;,, g;, Ableitungen 
nach x bis zur Ordnung 4, von wu die zweite Ableitung nach ¢ enthalt. Diese Ab- 
leitungen sind an geeigneten Zwischenstellen genommen. Wir setzen sie alle als 
existent und beschrankt voraus, im Sinne der in der Einleitung getroffenen Ver- 
einbarung tiber ,,hinreichend glatt’. Wir verzichten, den ausgerechneten lokalen 
Abbruchfehler hier wiederzugeben. 

Der Vollstandigkeit halber sei noch erwahnt, daB bei Auffassung von (42) 
als Beschreibung eines diskreten Irrfahrtprozesses mit Att folgende Bezie- 
hungen gelten: 


(Ax,r/At = a5 {0;4(¥ + 5 hyo t) —B;;(e— +h, #)} 


+4(> ae Sy ) + hy, t) — bj, (% —h,, t)} 
k=j+1 
=>) Pal) +09), 
j=l 


1 


1 ep Laat & 
Cx Anridt= 5 {0y4(*+ 4 byt) +0,)(x— 5 byt} — 4 (+ od 
{8jn(% +My, t) — 29,4 (x, t) +8;,(% —hy, O} 
= b,,(x) +0(h), 


(An; A %,>[At = (Bip (x + hej, t) + bj (6 — fej, #) +B jn(% + hp, t) +8 ;4(% — hy, O} 
=b,, (x,t) +O), 7+k. 


5. Wie notwendig sind die Diagonaldominanzbedingungen? 


Wir fragen, ob fiir beliebig kleines 4 Differenzenapproximationen nichtnega- 
tiven Typs mit dem Nachbarschaftssystem 9t(«; h) fiir die parabolische Differen- 
tialgleichung 
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existieren kénnen, wenn die abgeschwachte Diagonaldominanzbedingung 


j-1 n 
(43) b(n) =(D + > Jl jal A, a Oe Se 
k=1 k=j+1 
in einer einen Punkt (&, s) enthaltenden offenen Menge des IR" x {t]/OSt<T} 
verletzt ist, d.h., wenn dort fiir ein 7 =7 in (43) statt = das Zeichen < steht. 
Da man die Koordinaten x; umnumerieren kann, kénnen wir ohne Beschrankung 
der Allgemeinheit 7) =1 annehmen. Es geniigt, diese Frage fiir (8) zu untersuchen, 
da sowohl (FP) als auch (FPS) in diese Gestalt umgeformt werden kénnen. Ferner 
ist es geboten, das in (25) und (40) auftretende 6 durch 0 zu ersetzen, da in Sonder- 
fallen d=0 erlaubt ist (man vgl. den Abschnitt ,,Sonderfalle“ in Kapitel 3). 
Wegen der groBen Anzahl der Konsistenzgleichungen, zu denen die Nicht- 
negativitatsbedingungen noch hinzukommen, scheint das Problem fiir m =3 sehr 
schwierig zu sein. Der Fall 1 = 2 14Bt sich behandeln durch Modifikation einer von 
Greenspan fiir elliptische Differentialgleichungen entworfenen Methode ([8], 
S. 62—70). Es gilt der 


Satz 4. Die Differentialgleichung (8) mit n =2 set gleichmaBig parabolisch. Die 
b;,=5,;, 4; und € seven hinrerchend glatt®, speztell seien ste heschrankt. In einer 
einem Punkt (&, s) enthaltenden offenen Menge M sei 


(44) by1(%, t)<|d12(x, 2)|. 


Dann existiert in diesem Punkt fiir geniigend kleines h keine Differenzenapproxt- 
mation nichtnegativen Typs zu (8) mit dem Nachbarschaftssystem N(x; h) = 
{x + 9h, + Oghe|0.= —1,0, 1; 02= —1, 0, 1}. 

Beweis. Fir geniigend kleines # kdnnen wir annehmen, daB die konvexe Hiille 
von (é, h) x {t|s —tr<t<s-+t} in M liegt. Wir greifen den Punkt (€,s) aus 
M heraus, machen x = zum Zentrum einer Familie 3t (x; 4) von Nachbarschafts- 
systemen und beachten die Konstanz von o in (9). Eine konsistente Approxima- 
tion (10) nichtnegativen Typs muB die Bedingungen (K1), (K2), (K3), (NN) er- 
fiillen mit R= {7| «+ rheN(x; h)}. 

Wegen (K1) und (NN) sind die #, beschrankt, fiir eine geeignete Nullfolge {h,}, 
%=1,2,..., > 00, streben also die #,(x, s; h,, t) gegen Grenzwerte B,=8,(%, s). 
Dabei gehen (K1), (K2), (K3) iiber in 


(45) Prat, © Dir Bay D170 = 208, ,: 
reR reR rER 
(NN) geht tiber in 
(46) 6, Qindiny FEN. 
Setzt man 
(47) epee Boy =A, ate BL. fin) ore (05.0), 


wobei der einfache Index m der Reihe nach die Werte 1, 2, 3, 4,5, 6, 7,8 an- 
nimmt, wenn der Doppelindex 7’ der Reihe nach die Paare (1, 0), (0, 1), (—14, 0), 
(0, —1), (4,1), (—1, 1), (—1, —1), (4, —1) annimmt, so geht (45) tiber in das 


9 So glatt, wie es fiir die Erfordernisse der Konsistenztheorie erforderlich ist. 
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System 
+ Oy + Oy + Oy + Og + Og + Ks + Hy + KH, + XH =O 
apes — Os + &; —U% — OU, + Ug = O 
+ XH — 4, +5 + %% — h, —% =O 
(48) 


+ oy + Og + Os + Hg + Oy + Hy = 20 dy4 
+ a5 —% + 4, —g = 200, 
+ Oe + Oty + Xs + %g + a7 + %g = 20 dogo, 
das mit a, %7, % als Parametern eindeutig nach den GrdBen ap, a1, %, %y, %q, He 


aufgelést werden kann. Wir verzichten darauf, die vollstandige Losung anzu- 
schreiben, sondern begniigen uns mit 


(49) Oty = 00,1 — 20 by —%y +H, — 20g 
und 
(50) Og = 0D) —%q — OH. 


(49), (50) und die Ungleichungen 
(54) Ogee Op Oe EUR ee ea ce ree ce 


miissen erfiillt sein, wenn die Approximation von nichtnegativem Typ sein soll, 
was wir als Antithese zur Behauptung des Satzes 4 annehmen wollen. 

Ohne Beschraénkung der Allgemeinheit konnen wir nun voraussetzen, daB 
by (x, s)>0O ist. Denn wenn 0,.<0 ist, kann man (8) durch die Variablensub- 
stitution €,=—.x,, ,=%, in eine parabolische Differentialgleichung gleicher 
Struktur mit £;, anstelle von );, transformieren, ftir die aber im interessierenden 
Punkt B,.> 0 ist. 

Wir nehmen also an, es sei 0,. >0,, im Punkte (x, s). Aus (54) und (49) folgt 
dann 2%—a, +, S0 b,, —20b,,.< —ob,,, also o%>0b,,+2a,+a,, mithin 
wegen (51) 

(52) helio EtO aes 


Andererseits ist wegen (54) und (50) a, So ),, —a% Xo },,, also im Widerspruch 
zu (52) 

(53) A So by;. 

Damit ist der Satz bewiesen. 


6. Beispiele 


Die Satze 2 und 3 erfordern fiir 7 +-k die Angabe von Korrekturfunktionen g;,. 
Dies ist sehr einfach, falls keine der Funktionen 0;, ihr Vorzeichen wechselt. Falls 
die Funktion 6;, ihr Vorzeichen wechselt, bestimme man g;, in der Umgebung 
des Nulldurchgangs von b,, durch geniigend glatte Abrundung der ,,Kante“ des 
Graphen der Funktion | ;,|, dabei ist auf die Bedingung (24) bzw. (41) zu achten. 
Falls die Diagonaldominanz der Matrix B=(b;,) nur schwach ausgepragt ist, 
kann dies eine unangenehme Aufgabe sein. 
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Zur Illustration diene die Matrix (mit  =3) 


5 sin(x%,-+ %3) cost 
B(x, t) = | sin(x, + x.) $ e~* 
cost ase 4 
Aus (25) folgt 
0 Sn (x, t) 
< min {5 —|sin (x, + x)| —|cos#|, $ —|sin (x, + %5)| —e~', 4 —|cosé| —e~*. 


Eine mégliche Wahl ist 7 (x, t) = $=6. Diese Wahl geniigt auch der Bedingung 
(40). Fir die g;, betrachten wir nur die starkere Bedingung (24), die hier so lautet 


O54 (%, | Seje(* 2) S|0;.(% +3, y=F- 


Kritisch sind die Nulldurchgange von 6, .= sin (*,-+ %,) und b,,= cos ¢. Man kommt 
zum Ziel mit 


Bye (st) == 3 sin* (xg) gag (as Se, orig = FS cost te. 
Mit s=+*,+ x, ist nadmlich 
12 —| 2, =$+sin?s —|sins| = } —|sins| (1 —|sins}), 
Se3 —|bo3| = $ + cos?s —|coss| = } —|coss| (1 —|coss|) 
und fiir OSyX1 istOSy(1—y)S 
2S6ie—|he) Se, [he) +t Sor2S| 2] +2, 


(24) ist also erfiillt fiir g,,. Analog zeigt man, daB (24) fiir g,, erfiillt ist. 


<}, mit y=|sin s| folgt also 


AbschlieBend sei fiir den Fall » =2 ein numerisches Beispiel betrachtet, und 
zwar eine inhomogene Gleichung, deren Lésung bekannt ist: 


j=1 k=1 
mit 
ie x 3 sin (x, + %9) 
sin (x, + %») 2 ; 


= 2¢{2 + sin (x, + %)} + (2) (1 +2) {7sin (x, + %2) — 4c0s (2 (x, + %9))} 
und der Losung 
u(x, t) = (1+#) {2+sin(*,+ x,)}. 


Die Gleichung hat selbstadjungierten elliptischen Anteil, und wir diskretisieren 
gemaB Kapitel 4 mit der Modifikation (10) (mit /=/). Wir nehmen 7 (x, t) = 
3 —|sin(x,+ %,)| und, vertraglich mit (41), g2(%,#)=1. Damit finden wir mit 
den Formeln (32), (34), (37) 

p(%,%;t)=1—S50, p(x, x+hst)h=o/2, p(%,x the; t)=o, 
p(x, x + Ah, + phy; t) = (0/4) {Aw (sin (x, + x2 + AA) +sin (x, + % +h) +2}, 
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und haben, mit o<} und *=mh, 


U(x, 0) = g(x) =2+sin (x, + %), 
U(x,t-+7)= >) p(z, x32) U(z,t) +77 (2,2). 
) 


ZEN (x 

Bei der Rechnung beachten wir die Periodizitat der Diffusionsmatrix B (x, t) 
und der Anfangsbedingung und beschrénken uns auf das Quadrat —mS% Sa, 
— < %, <x mit periodischer Fortsetzung der p (z, x; ¢) und der jeweils erhaltenen 
Werte U(x, 2). 

Die Rechnung!® wurde durchgefiihrt auf der IBM 360/91 des Instituts fiir 
Plasmaphysik mit ungefahr 16 signifikanten Dezimalstellen, und zwar mit den 
Schrittweiten h=2/16 und h=z/32 und dem Parameter o =0,15. Fiir die Werte 
t=2 und ¢=4 und die Stellen x= (j,7/4, 7.7/4), 71, 72 gamzzahlig, wurden die 
Werte U —u fiir beide Schrittweiten ermittelt und durcheinander dividiert. Da 
der Fehler asymptotisch wie O(h?) abnimmt, miissen diese Quotienten in der 
Nahe von 4 liegen. Infolge der Symmetrie-Eigenschaften der speziell gewahlten 
Differentialgleichung und ihrer Lésung treten in jedem Resultat-Schema nur fiinf 
voneinander verschiedene Eintrage auf, die beispielsweise an den Stellen x= 
(0, 7222/4) stehen, wobei 7, = —2, —1,0, 1,2. Die gerundeten Werte geben wir 
in Form einer Tabelle: 


t 45 Umith=n2/16 Umith=a/32 u Fehlerquotient 
2 24 5,01079 5,00269 5,00000 4,01 
2 4) 6,46350 6,46427 6,46447 4,84 
2 0 10,02634 10,00668 10,00000 3,94 
2 1 13,62725 13,55841 135538553 4.01 
2 2 15,12501 15,03109 15,00000 4,02 
4 2 17,00923 17,00226 17,00000 4,08 
4: =—4 21,93758 21,96889 21,97918 4,04 
4 O 34,02561 34,00671 34,00000 3,81 
4 4 46,24131 46,07578 46,02082 4,01 
4 2 51,32175 51,07997 51,00000 4,02 


Der in der zweiten Zeile auftretende ungewoéhnlich groBe Fehlerquotient 4,84 
hat offenbar die Ursache, daB U —wu dort ein anderes Vorzeichen hat als an den 
anderen Stellen (x, ¢), fiir die Resultate angegeben sind. 


Zusammenfassung 


Fiir die partielle Differentialgleichung von Fokker und Planck im IR” (Vor- 
warts-Differentialgleichung von Kolmogorov) mit orts- und zeitabhangiger Drift 
und Diffusionsmatrix werden, unter einigen zusatzlichen Bedingungen, explizite 
Differenzenschemata angegeben, die zwei wesentliche Eigenschaften des beschrie- 
benen Diffusionsprozesses imitieren, ndmlich Nichtnegativitatserhaltung und Sub- 
stanzerhaltung. Diese Schemata kénnen auch interpretiert werden als Beschrei- 
bung diskreter instationérer inhomogener Irrfahrten. Sie sind stabil in der 


10 Fiir die Programmierung bin ich Herrn J. Steuerwald zu Dank verpflichtet. 
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Maximumnorm, und man kann mit ihrer Hilfe bei gegebener Anfangsbedingung 
die Differentialgleichung numerisch approximativ lésen oder aber den durch sie 
beschriebenen Diffusionsproze8 approximativ simulieren (Monte-Carlo) und ver- 
anschaulichen. Wesentliche Bedingungen sind die Beschranktheit der Koeffi- 
zienten der Differentialgleichung und einige ihrer Ableitungen und eine starke 
Diagonaldominanz der Diffusionsmatrix. Die Resultate lassen sich verallgemei- 
nern auf allgemeine lineare parabolische Differentialgleichungen (Substanzerhal- 
tung ist dann i.allg. nicht mehr vorhanden). Die beschriebene Methode erlaubt 
auch, fiir lineare elliptische Operatoren im IR” Approximationen nichtnegativen 
Typs anzugeben, wenn diese analogen Bedingungen geniigen. 
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Eine Méglichkeit zur Konvergenzbeschleunigung 
bei Iterationsverfahren fiir bestimmte nichtlineare Probleme 


DIETRICH BRAESS 


Eingegangen am 26. Februar 1969 


Nichtlineare Gleichungssysteme, nichtlineare Approximationsaufgaben und andere 
nichtlineare Probleme behandelt man oft numerisch mit Iterationsverfahren. In jedem 
Iterationsschritt ist eine lineare Aufgabe zu lésen, die eine Naherung des vorliegenden 
Problems darstellt. Solche Verfahren lassen sich in einem allgemeinen Rahmen formu- 
lieren. Eine Méglichkeit zur Konvergenzbeschleunigung bietet sich dann bei Problemen, 
die in einem bestimmten Sinne teilweise linear sind, indem der klassische Algorithmus 
durch einen kleinen OptimierungsprozeB erweitert wird. Fiir Probleme mit Exponen- 
tialsummen erhalt man diese Form der Beschleunigung auch mit Hilfe von Invarianz- 
betrachtungen. 


1 


Nichtlineare Probleme behandelt man oft mit Hilfe von Iterationsverfahren. 
Mit jedem Schritt wird eine lineare Aufgabe numerisch gelost, die als Naherung 
des zugrunde liegenden Problems verstanden werden kann. Das bekannteste Ver- 
fahren von diesem Typ ist das Newtonsche Verfahren zur Losung nichtlinearer 
Gleichungssysteme. 


F(@j]E (my, %..%)=0 Ce yok |S (1) 


Die Gleichungen fassen wir durch Einfiihrung der -wertigen Funktion F [«] zu- 
sammen. 


EF Vee) =O. 


Zu der Naherung beim v-ten Iterationsschritt «” wird die (differenzierbare) Funk- 
tion F [«] durch die lineare Funktion 


F [a] =F [a”] + (« —0)' VF (a?) (2) 


formal ersetzt, und als nachste Naherung «”** wird der Lésungsvektor des linearen 
Systems 


F(a) Ole = 41, 2.2) (3) 


betrachtet. Nach dem gleichen Prinzip gelangt man zu Newtonschen Verfahren 
fiir andere Aufgaben, z.B. fiir die Anpassung nichtlinearer Funktionen nach der 
Methode der kleinsten Quadrate [1,6]. Das Konstruktionsprinzip formulieren 
wir im nachsten Abschnitt so allgemein, daB sich u.a. auch die Remes-Algorithmen 
fiir die nichtlineare Tschebyscheff-Approximation [7, 10, 11] einordnen. 
Bekanntlich braucht bei solchen Verfahren in vielen Fallen keine Konvergenz 
einzutreten, wenn die Ausgangsnaherung nicht schon hinreichend nahe bei der 
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Lésung gewahlt ist. Eine Méglichkeit zur Konvergenzbeschleunigung bietet sich 
bei ,,teilweise linearen‘‘ Problemen, d.h. wenn die Funktion F [a] von einigen 
Parametern linear abhangt. So sind z.B. bei der Exponentialapproximation die 
Funktionen 
l 
F(%, ay =" ee * (7 = 21) (4) 
k=1 
linear in den ersten / Parametern. Bei der Anwendung der allgemeinen Vorschrift, 
die im dritten Abschnitt beschrieben wird, sind bei der modifizierten Iteration die 
linearen Parameter? durch einen OptimierungsprozeB zu bestimmen, und nur die 
iibrigen Parameter? werden in der tiblichen Weise berechnet. Im vierten Abschnitt 
folgen Invarianzbetrachtungen fiir die Iteration bei der Exponentialapproxima- 
tion. Im fiinften Abschnitt werden noch zwei weitere Beispiele fiir die Anwendung 
des Konstruktionsprinzips gegeben. 
Der Erfolg der modifizierten Iteration zeigte sich gerade bei der Anpassung mit 
Exponentialsummen sehr deutlich, weil das iibliche Newtonsche Verfahren bei 
mehrgliedrigen Summen wegen der Rundungsfehler meistens versagt [2]. 


2 
Ehe wir die Verfahren allgemein formulieren, werden die Voraussetzungen an 
Hand eines Algorithmus fiir die nichtlineare Tschebyscheff-A pproximation erlau- 
tert. Gegeben sei eine Menge von Funktionen F [a] =/’*(x, «), die auf einem reellen 
Intervall X definiert und die durch den Parametervektor « charakterisiert sind. 
Zu einer vorgelegten Funktionen /(x) ist fiir den Anpassungsfehler 


P(x) =|F [a] —f|:= max | F(x, «) —/(x)| (5) 


das Minimum zu ermitteln. Die Bestimmung des optimalen Parametervektors ist 
nun selbst dann nicht in einem Schritt mdglich, wenn F[«] linear von den Para- 
metern abhangt. 

Eine starke Vereinfachung ergibt sich, wenn man die Anpassung nur tiber 
einer endlichen und diskreten Punktmenge {x;} und nicht iiber dem ganzen Inter- 
vall ausfiihrt. Im linearen Fall ist die Minimierung von 

max |F(x;, «) —f(%,)| (6) 
fiir endliche Punktmengen {x,;} < X mit Hilfe eines linearen Programms méglich, 
das sich in Spezialfallen sogar auf ein lineares Gleichungssystem reduziert. 

Bei der Konstruktion der besten Approximation nach Remes wahlt man ab- 
hangig von der Naherung «” eine Menge X” c X mit (mindestens) » +1 Punkten, 
die im allgemeinen Extremalpunkte der Fehlerkurve F(x, «’) —/(x) sind. Die 
Fehlerfunktion 

@" (x) = max | F(x, «) —f(x)| (7) 


xEexv 


besitzt, wenn man sie auch als Funktional von F auffaBt, folgende Eigenschaften. 


1 Zum Beispiel sind das in (4) die Parameter a ... a. 
2 Zum Beispiel sind das in (4) die Parameter «),,... %;. 
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I, Fir a= a gilt 


II. Zu jeder linearen Funktion 
I 
F (a]=Hy+ 2 % 7, 
k=1 


1aBt sich numerisch das Minimum von ®’(«) bestimmen. Dabei ist fiir die Anzahl 
der (freien) Parameter «, jede Zahl / <n zulassig. 


Unter diesen Voraussetzungen la4Bt sich die Iterationsvorschrift fiir Fréchet- 
differenzierbare Funktionen F [«] allgemein formulieren. 
4. Fiir «=«” werden zum Aufbau der ,,linearisierten Funktion“ F die Ablei- 
a : , 
tungen a, 2 {a”] ermittelt. 


2. Wenn in dem Funktional ©” an Stelle von F die Naherung F eingesetzt 


wird, erhalt man die Funktion PD’ (x). Fiir diese ist das Minimum « =& zu ermit- 


teln. Es sei wi RS (9a) 


die berechnete Parameteranderung und 
A® := D(a’) — B (&) = D(a’) —B’ (4) (9b) 
die fiir das linearisierte Problem erreichte Verbesserung. 
3. Es werde zundchst der Dampfungsfaktor ¢=1 gesetzt und getestet, ob 


D(e’ +¢tAa) S D(a’) —4tA®@ (10) 
gilt. Sofern die Bedingung (10) erfiillt ist, wird ein neuer Iterationsschritt mit 
oe tte’ +tAa (411) 


begonnen. Andernfalls wird der Dampfungsfaktor ¢ halbiert und der Teilschritt 3 
wiederholt’. 


Die Iteration erzeugt eine monotone Folge 
P (0°) = O(ot) > Oe) =S---. (12) 


Sie wird beendet®, sofern wegen der Rundungsfehler keine Verbesserung mehr zu 
erreichen ist. 

Bei den klassischen (Newtonschen) Iterationsverfahren gibt es nur die ersten 
beiden Teilschritte, und es wird direkt «’**=a%=«" +A gesetzt. Das entspricht 
einem Dampfungsfaktor ¢= 14. 

Die Monotonie (12) wird durch den Test im dritten Teilschritt erreicht. An 
die Stelle des Tests tritt bei einigen Autoren die Minimierung von ®(«’ +¢ 4a) 
zur Bestimmung des Faktors ¢ [4, 6]. Bei einer anderen Variante wird der Faktor 

3 Da8 in Programmen ftir Rechenmaschinen vorsichtshalber fiir die Zahl der 


Iterationsschritte und Teilschrittwiederholungen eine Schranke einzubauen ist, lassen 
wir in dem formalen Programm auBer acht. 


4 Die im folgenden dargestellte Modifikation zur Konvergenzbeschleunigung kann 


auch ohne Dampfung (d.h. fiir =1) verwandt werden, wenn auch dort der Test mit 
seinen IXonsequenzen unterdriickt wird. 
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mit Hilfe der zweiten Ableitungen so abgeschitzt, daB eine zu (10) aquivalente 
Relation erfiillt ist [8]. In jedem Fall ist die Iteration in elektronischen Rechnern 
wegen (12) leicht zu kontrollieren, da im Gegensatz zu den klassischen Verfahren 
die bekannte Kafigbildung und andere Oscillationen [9] nicht auftreten. Die 
Konvergenz erfolgt in einer Form, wie sie von Ostrowski fiir Gradientenverfahren 
beschrieben wurde [8]. 

Eine Gefahr dafiir, daB die Iteration vorzeitig abgebrochen wird, besteht 
allerdings bei kleinen Dampfungsparametern. Denn @(«’) und ®(«’ +¢ 4) unter- 
scheiden sich dann nur wenig, und die gerundeten numerischen Werte erfiillen 
nicht mehr die Relation (10). Aus diesem Grunde erwies sich eine Modifikation 
zur Konvergenzbeschleunigung als zweckmaBig. 


3 


Das im letzten Abschnitt beschriebene Konstruktionsprinzip modifizieren wir 
bei Problemen, bei denen die Funktion F{«] von einigen Parametern «, linear 
abhangt, also in der Form (13) darstellbar ist 


1 
Fla]= )'%,H, [Xp Hye +++ %y] (<n). (13) 


k=1 


Um diese Struktur noch deutlicher hervorzuheben, teilen wir die Parameter in 
zwei Gruppen und schreiben mit «= (A, y) 


l 
Fla]=F[B, y] = Bi [y] 
und entsprechend 


D(x) =D(B,y), D(a) = PB, y). 


Die Berechnung neuer Naherungswerte fiir die Parameter wird in den zwei Grup- 
pen unterschiedlich vollzogen. Bei der ,,modifizierten Iteration“ ersetzen wir den 
Teilschritt 3 aus der ,,einfachen Iteration‘ durch die Vorschrift: 


3’. Es werden «&” und Aa gemaB «” = (f’, y”) bzw. Aa = (AB, Ay) aufgespalten, 


und es wird t= 1 gesetzt. Fiir f werden die Werte f ermittelt, bei denen 
D"(B, y" +t Ay) 


minimiert wird, und es wird getestet, ob 
O(, y" +tAy) < Ow) —114@ (10’) 


gilt. Sofern die Bedingung (10’) erfiillt ist, wird ein neuer Iterationsschritt mit 
oe (B, y’ +t Ay) begonnen. Andernfalls wird der Faktor ¢ halbiert und der 
Teilschritt 3’ wiederholt. 

Bei der erweiterten Vorschrift wird nicht nur mit jedem Jtevationsschritt ein 
Minimumproblem geldst, sondern auch mit jedem Démpfungsschritt. Allerdings 
ist die Anzahl der freien Parameter bei den Dampfungsschritten reduziert, und 
der Arbeitsaufwand ist dementsprechend kleiner. Andererseits gilt offensichtlich 


® (B, y’ +tAy) < © (B’ +t AB, y’ +t Ay), 
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und so ist fiir manche Werte von ¢ wohl die Bedingung (10’) aber nicht (10) erfiillt. 
Fiir die Iteration nach der urspriinglich gegebenen Vorschrift waren deshalb noch 
weitere Dampfungsschritte erforderlich, ehe eine bessere Naherungslésung gefun- 
den ist. Damit wird der Mehraufwand bei der Ausfiihrung von 3’ gerechtfertigt. 


4 

Es scheint zundchst so, als ob durch die Aufteilung der Parameter in zwei 

Gruppen und durch die verschiedene Behandlung eine Symmetrie der klassischen 

Verfahren zerstért wird. Aber es zeigt sich im Gegenteil, daB bei der Anpassung 
mit Exponentialsummen 


1 
=> Pom (14) 


gerade durch die Modifikation eine Symmetrie zwischen den Parametern f, und 
y, erreicht wird. In Hinblick darauf diskutieren wir das Verhalten der Iterations- 
folgen bei Verschiebung des Approximationsintervalls. 

Uber einem reellen Intervall X sei die stetige Funktion f(x) durch Exponential- 
summen (14) im Tschebyscheffschen Sinn bestméglich zu approximieren. Sei & 
eine feste reelle Zahl und 

X = {&|%=x4+6, x€X} 


das transformierte Intervall. Wir fiihren einen Verschiebungsoperator 


PC (=e) 
durch (15) 
(Pf) (x) =f(% —€) 


ein. P bildet die Familie der Exponentialsummen in sich ab. Offensichtlich gilt: 
Wenn F [«] beste Approximation zu / im Intervall X ist, dann ist auch PF [«] 
beste Approximation zu P/ im Intervall x 

Definition. Sei F [«”] eine Iterationsfolge zur Konstruktion der besten Appro- 
ximation zu f/€C(X) bei der Anpassung in X, und F [&”] sei die Folge zu Pf bei 
der Anpassung in X. Die Iteration heiBt invariant gegeniiber Translationen, wenn 
aus 

Ele} = PF ile) 
die Relationen 
FP late PR iat firin 12 53he. 

folgen. 


Sel num 2'(4;-0°). == 2 B, e’”** die Ausgangsnaherung und ¢(x) =/(x) —F (x,«°) 


die Fehlerfunktion. Teen (2) und (6) bedeutet die Minimierung von ®°(«) die 
Bestimmung von Anderungen J «,, so daB an den Punkten x, fiir die Differenz 


— >) aes Fae Fw, a) 


das Maximum der Betrage méglichst klein wird. Wie man durch Ausrechnen der 
Ableitungen fiir die Exponentialsummen (14) erkennt, ist dazu die Funktion e°(x) 
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der Form 


mit den freuen Parametern o;, t, aufzusuchen, welche die Fehlerkurve ¢(x) am 
besten approximiert. Aus e°(x) ergeben sich die Parameteranderungen 


AB,=o,, Ay,=T/Bs (B, +0), R=1, 2,...,0. (17) 


Fir das transformierte Problem sei (16) erfiillt. Dann hat 
Lae l 
F [6°] = PF [o®] = > B, e*** aa (B,er er (18) 
k=1 k=1 


die gleichen Frequenzen y, wie E [«°] und die zugehérige Fehlerkurve ist ¢ [&°] = 
Pe{«®}. Wir kénnen annehmen, daB beim linearisierten Problem die Anpassung 
auf den um die GréBe € verschobenen Punkten 


~— x; téeX 


geschieht; denn die Wahl hangt bei Remes-Algorithmen [7] nur von der Gestalt 
der Fehlerkurve ¢(x) ab. Also lautet die Losung des zweiten Teilschrittes 


1 
é*(«) = (Pe) ( =k (G, +, x) evF* 

mit 7 
6, NO ere en tS ea (19) 
Aus (17), (48) und (19) folgt unmittelbar, daB bei beiden Iterationen die Frequenz- 
anderungen tibereinstimmen, daB sich aber die Anderungen der Faktoren f, nicht 
entsprechen. Die einfache Iteration fiihrt also nicht auf invariante Folgen. Bei 
der modifizierten Iteration werden dagegen zu den vorgegebenen Frequenzen 


jeweils die optimalen Faktoren #, eingesetzt. Da dieser Vorgang invariant gegen- 
iiber den Verschiebungen ist, gilt das auch fiir die ganze Iterationsfolge. 


5 

Es seien noch zwei Beispiele fiir Iterationsverfahren genannt, bei denen das 
beschriebene Konstruktionsprinzip angewandt wird. Es geniigt, fiir die Beispiele 
das Funktional ® zu definieren und fiir den 2. Teilschritt die spezielle Vorschrift 
zur Berechnung von 4a und A® anzugeben. Die anderen Teilschritte sind schon 
durch den allgemeinen Formalismus vollstandig festgelegt. 

Beispiel 1. Anpassung nach der Methode der kleinsten Quadrate [1, 4, 7]. Es 
sind fiir den Funktionsansatz F(x, «) die Parameter «, so zu bestimmen, daB die 
gewichtete Fehlerquadratsumme 


= > 0; [F(%;,«) —}%, w,>0,m2n (20) 
j=1 
minimal wird. Es wird ®’(«) = ®(«) gesetzt und A« ergibt sich aus dem Glei- 


chungssystem 
JWjAa—J'Wh=0. 
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Dabei berechnet sich die Jakobische Matrix J, die Wichtematrix W und der 
Fehlervektor h gemaB 


p 
Jin = Gag E(t ©), 


k 
h,=;—F'(x;, &). 


Ferner gilt AD=Aa' J'Wh. 
Beispiel 2. Lésung nichtlinearer Gleichungssysteme [3, 5]. Eine Bewertung 
von N&aherungslésungen kann man wie in [5] durch die Fehlerquadratsumme 


oder durch das Fehlermaximum 


® (a) = max | F,(2)| 
einfiihren. Dann ist auch ®’(«) = ®(«) modglich. In beiden Fallen fiihrt die Mini- 
mierung im zweiten Teilschritt auf das zu (3) aquivalente Gleichungssystem 
VF [a] -da= —F [a]. Wegen ®’ («” +4) =0 vereinfacht sich die Forderung (10) 
zu 


Pa? +140) S$ (1-3) Oe), 
und eine entsprechende Vereinfachung ergibt sich fiir (10’). 


6 


Zum Schlu8 betrachten wir noch ein numerisches Beispiel. Die Untersuchung 
der Reaktionsgeschwindigkeit bei bestimmten chemischen Prozessen fiihrt auf die 
Differentialgleichung 


IIE ay 
du au—by+c’ 


u>dO. 


Die Parameter a, b und c sind aus den MeBwerten zu bestimmen. Die Lésungen 
der Differentialgleichung haben die Gestalt 


i + + sey & MER -é* mit «=e7*. 
Dies legt den Ansatz 
y=a, +a,¢~*+a,6%* 


nahe. Jedoch wird dann die Jakobische Matrix fiir «,—+0O*singular. Deshalb ist 


der aquivalente Ansatz 
ek *% — 4 


Y= +0, *--O, 


oq 


geeigneter. Der Einflu8 der Beschleunigung zeigt sich nun deutlich bei der Mini- 
mierung der Fehlerquadratsumme (20) fiir die Punkte (Testwerte) 


4, = 1 VV as ¢=1,2...7 
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mit den Gewichten W,=1. Als optimale Anpassung erhalten wir die Funktion 


— 0,0867 4 __ 4 


v=1,91 —1,10e-* +0,763 S568) S 


Die herkémmliche Iteration verlangt 30 Schritte. Wenn man aber ausnutzt, daB 
die Funktion in dem ersten (bzw. in den ersten beiden, bzw. in den ersten drei) 
Argumenten linear ist, reduziert sich die Schrittzahl auf 21 (bzw. 18, bzw. 12). 


40. 


ale 


Herrn Professor Dr. H. Werner méchte ich fiir seine Anregungen danken. 
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Existenz und Eindeutigkeit 
fiir Losungen nichtlinearer Randwertprobleme 


WALTER PETRY 


Eingegangen am 19. Juni 1968 


1. Einleitung 
Ehrmann [3] und Conti [1] geben Existenzsatze fiir die Losung des Problems 


(4.1) L,u=M,u, 
(1.1 b) L,u=M,u 


an, wobei L, und L, lineare Abbildungen eines Banachraumes B, in Banachraume 
B, bzw. B, und M, und M, (nichtlineare) Abbildungen von By in B, bzw. By sind. 
Hierbei wird vorausgesetzt, daB die Operatoren M, und M, auf dem ganzen Raum 
B, beschrankt sind, d.h. es gibt positive Konstanten a, a,, so daB fiir alle we By 
gilt 

(1.2) |M, ul, Sa, ||Mzullp, Say. * 

In einer spateren Arbeit von Ehrmann [4] wird die Beschranktheitsbedingung (1.2) 
durch eine Bedingung der Form 


(1.3) Mims, Sa +, (ele,  [Mamlle, S42 +52 lle, 


ersetzt, wobei 0; (¢=1, 2) nichtnegative Konstanten sind, die gewissen Bedin- 
gungen geniigen miissen. 

In dieser Arbeit werden unter der Voraussetzung, daB der durch Lju=O, , 
L,u=O,. erzeugte lineare Operator formal selbstadjungiert ist, zwei Existenz- 
satze und zwei Existenz- und Eindeutigkeitssatze fiir ,,verallgemeinerte“’ Lésun- 
gen des Problems (1.1) bewiesen, welche die lineare Beschranktheit von M, nicht 
benotigen. 

Diese Satze ergeben sich aus einem abstrakten Existenzsatz und einem ab- 
strakten Existenz- und Eindeutigkeitssatz des Verfassers [8]. 

Als Anwendung der erhaltenen Ergebnisse betrachten wir fiir eine offene be- 
schrankte Menge 2c R? die Differentialgleichung 


(1.4) —Au(s)=f(s,u(s)) fir seQ 


mit gewissen nichtlinearen Randbedingungen fiir «(s) und mit auf $ = {(s,)|s€Q, 
u€R}} stetigem /(s, w). Ebenso untersuchen wir die Differentialgleichung 
(1.5) —u''(s)=g(s,u(s)) fir se]o,4[ 


1 ||: \|z bezeichnet die Norm von B. 
2 Op bezeichnet das Nullelement von B. 
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mit gewissen nichtlinearen Randbedingungen, wobei g(s, #) am Rand Singulari- 
taten besitzen darf. Fiir s€]0, 1[ sei g(s, ) stetig. 
Die Satze von Ehrmann und Conti sind auf diese Beispiele nicht anwendbar. 


2. Hilfsmittel 


In diesem Abschnitt formulieren wir zwei abstrakte Satze, die in einer friiheren 
Arbeit des Verfassers [8] bewiesen wurden, und welche beim Beweis der Existenz- 
und Eindeutigkeitssatze fiir das Problem (1.1) bendtigt werden. 


Bezetchnung 1. Im folgenden bezeichne H einen Hilbertraum mit dem Skalar- 
produkt (-, -) und Z einen Banachraum. Ferner sei fiir 0Sa<co K,={y|yeB, 
lyle So} und Ly = » la Se}. 

Hilfssatz 1. Es seien die folgenden Voraussetzungen erfiillt. 

) Es existiere ein R mit 0S R< «~, so daB gelte: 


a) I(-,-) sei eine stetige Abbildung von H x Kp in B. 
Bp) Fir ye Kp sei S(-, y) eine stetige Abbildung von A in sich. 
ely Furs 6H sei.S (x, ) eine stetige Abbildung von Kp in H. 


B) Fiir yeK, und %,, *,€H gelte 
Re (S(x,y) —S(%2, 9), %1— 4%) 2C | % — Xelli 


(A 
( 
( 
( 
( 
) 


(2.4 
mit einer Konstanten C >0. 
(C) Fir yeK, gelte 
SS |S, le SM 
mit einer Konstanten M > 0. 
(D) Fiir x€Ly,c und yeKe gelte 
(2.3) (T(x, y[pSR. 
(E) T(-,-) sei eine kompakte Abbildung. 
Dann existiert mindestens ein x*, y* mit «*€Ly,c und y*€Ke von 
(2.4) S(%y)=On, y=T(x,9). 
Bemerkung 1. Setzt man S(x, y) = «x —Sp(x, y), dann kann das Problem (2.4) 
in der Form 
(2.5) x= So(4%,y), y=T(x, 9) 
geschrieben werden. 
Hilfssatz 2. Es seien die folgenden Voraussetzungen erfiillt. 
(A’) Es existieren Konstanten R und M mit OS R< wo undOSM< ~~, so 


daB gelte: 
(x) T(-, -) sei eine stetige Abbildung von Lyjc X Kz in B. 
(6) Wie (If) von Hilfssatz 1. 
(y) Es existiere eine Konstante L > 0, so daB fiir x€Ly,c und 4, yee Kp gelte 


(2.6) [So (%, M4) —So(*, Ve ) ler SL | Vrs Valle - 
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(B’) Wie (B) von Hilfssatz 1. 
(C’) Wie (C) von Hilfssatz 1. 
(D’) Wie (D) von Hilfssatz 1. 
(E’) Fiir *,€Ly/c und y;eKp (¢=1, 2) gelte 


(2.7) |Z (41, 1) —T(%2, ve) lp Sfx la — %alle + Fe || 1 — Yelle 
mit Konstanten k; 20 (1 =1, 2) und 
16. 
(2.8) ri= hye the<i. 
Dann besitzt das System 
(2.9) %=So(%, ¥), y=T(x, y) 


genau eine Lésung x*, y* mit y*e€Kp. Es ist x*€Ly,~. Ferner konvergiert fir 
yo €K p das Iterationsverfahren 


(2.10) X41 =So(%pias My)» Veti=L(% 41%») (v= 0,1, 2-.-.) 
gegen die eindeutige Losung x*, y*. 


3. Existenz- und Eindeutigkeitssatze 


In diesem Abschnitt untersuchen wir Existenz und Eindeutigkeit einer ,,ver- 
allgemeinerten‘“‘ Lésung der nichtlinearen Randwertaufgabe 


(3.14) (L, 4) (s)= (4,4) (s) fiir se€Q, 
(3.1 b) L,w=M,u. 
Hierbei bezeichnet 2 eine offene, beschrankte, meBbare Menge aus Pa (J =1) mit 
dem MaB |Q]. 

Fir die anschlieBenden Betrachtungen werden folgende Raume eingefiihrt: 

E”: Raum der reellen »-dimensionalen Vektoren mit dem Skalarprodukt 
la IE qh: do und der Norm ||z||zn = (<z, z>)* (Hilbertraum). 

LPG: Sere der auf 2 reellen meBbaren Funktionen z(s) mit ed |z(s)|?ds< 00 
(4 <p < oo) und der Norm |z|,»= (Jf lz | 2( s)|Pas). Hierbei mentee Funktionen, 
die fast iiberall auf Q gleich sind, als identisch betrachtet. 

f°: Es sei p=(,,..., 2,), dann sei ZG? = I L?‘(Q) mit der Norm |u|,» = 
bs Jel.) ( (Banachraum). : 


Fir p;=2 (¢=1, ..., m) erhalt man speziell einen Hilbertraum, der im folgen- 
den mit #* bezeichnet wird. 
V oraussetzung 1 


(I) Es sei L, ein linearer Operator mit Definitionsbereich D(L,)¢¥? mit p; =2 
(1 =1,...,) und Wertebereich R(L,) aus einem linearen Funktionenraum B. 
M, sei ein (nichtlinearer) stetiger beschrankter Operator von ¥? in #4 mit q= 
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(Cree) tnd: ka + os =41 (¢=14,...,). LZ, sei ein linearer Operator mit 
D(L,)¢D(L.)¢cH? und R(L,) aus einem Banachraum B. M, sei ein (nichtlinearer) 
stetiger beschrankter Operator von ¥? in B. 


(II) Sei By eine nichtleere Teilmenge von #90B und sei fEB,, dann existiere 
eine Abbildung K, so daB aus 


(3.2a) u(s)=(Kf)(s)  (s€Q) 
folgt w€D(L,) und 
(3.2b) (L,u)(s)=f(s) (seQ), L,uw= Og. 


Die Abbildung K sei linear und beschrankt von #7 in Y?. AufgefaBt als Abbildung 
auf ¥? sei K selbstadjungiert und positiv definit. 


(III) Fiir z€¢Y" gelte die Darstellung 
(3.3) Raa OFF 2, 


Hierbei sei ¥ ein linearer beschrankter Operator von #? in #? und #* der zu# 
adjungierte Operator, d.h. #* ist linear und beschrankt von #1 in #2. 


(IV) Fir y¢B existiere u(y) €D(L,) mit 
(3.4) (L,%(y))(s)=0 fir seQ. 


Hierbei sei w#,) ein linearer beschrankter Operator von B in #?. 


(V) Es existiere auf B ein linearer beschrankter Operator £ mit 
(3.5) Eyu (Ly) =L(L a) y=, 
d.h. L=(L, u)+. 

(VI) Fir w,e2? (7 =1, 2) gelte 
(3-6) JM (wt) (s) — MG (uta) (s), 4 (s) — Ma (s)> ds Sau, —glgn. ° 
Hierbei sei 


(3.7) aZ=0, C:=1—al#| 


gas 

(VII) Es existiere eine Konstante b> 0, so daf fiir we Y? gelte 
(3.8) |], up SO. 

(VIII) Z M, sei kompakt von Y? in B. 


Voraussetzung 2 

(I’) Es sei L, ein linearer Operator mit D(L,)¢ Bc ¥? und R(L,) aus einem 
linearen Funktionenraum B. M, sei ein (nichtlinearer) stetiger und beschrankter 
Operator von B, in B,c ¥*. Ly sei ein linearer Operator mit D(L,)¢D(L,) cB, 
und R(L,) aus einem Banachraum B. M, sei ein (nichtlinearer) stetiger beschrank- 
ter Operator von B, in B. 


3 Es ist £Pc $2 wegen p;=2 und |Q|<ov. 
4 |aljg-,, bezeichnet die Operatornorm von « als Abbildung von f in y. 
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(II’) Sei B, eine nichtleere Teilmenge von BoB und sei f€ B,, dann existiere 
eine Abbildung K, so daB aus 


(3-9) u(s)=(Kf)(s) (se) 
folgt weD(L,) und 
(3.9) (L,u)(s)=f(s) (se), Law =Oz. 


Hierbei sei K ein linearer beschrankter Operator von B, in B,. Es existiere ein 
linearer, beschrankter, selbstadjungierter, positiv definiter Operator K auf # 
mit J, K=K J,, wobei J; (i =1, 2) die Injektionsabbildungen von B; in ¥? be- 
deuten. 

(III’) Die Operatoren K} und Kj seien linear und beschrankt von ¥Y? bzw. B, 
in B, baw. Y? und es gelte Kt = J, K} bzw. K}= K*J,, wobei K? den eindeutig 
bestimmten selbstadjungierten, positiv definiten Wurzeloperator von K bezeich- 
net. 

(IV’) Fiir y€B existiere u(y) €D(L,) mit 


(3.10) (L,%(y))(s)=0 fiir seQ. 
Hierbei sei #,) ein linearer beschrankter Operator von B in B,. 

(V’) Wie (V). 

(VI’) Fir u,€B, (7 =14, 2) gelte 
(3.11) (JoM, (4) —JoM, (us), Jr — J, Uz) Sa Tim —J, ual. 
Hierbei bezeichne (-, -) das Skalarprodukt in #?, und es sei 
(3.12) a=0, C:=1—a|K+# fy, ,~>O0. 

(VII’) Es existiere eine Konstante b> 0, so daB fiir w€ B, gelte 
(3.13) |M, #\|5 Sb. 

(VIII’) LM, sei kompakt von B, in B. 

Lemma 1. Es sei die Voraussetzung 1 erfiillt. Ferner sei v*€ ¥?, y*€ B Lésung 


von 


v=H* M, (u(y) +2), 


(3.14) . 

v= LM, (u(y) +H). 
Dann gilt 
(3.15) u* := Uy (y*) +#v* Eg? 


und u*, y* geniigt den Operatorgleichungen 
u = Uy (¥) +4 M, (u) 


(3.16) 
Ly u(y) = Mu. 


Beweis. Die Behauptung folgt aus (3.14) unter Beachtung von (IV), (I), (II) 
und (V). | 
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Lemma 2. Es sei die Voraussetzung 2 erfiillt. Ferner sei v*€ #2, y*€ B Lésung 
von 


v= K}M, (u(y) + Kyu), 


(3.17) é : 
y=LM, (29 (¥) + Kv). 

Dann gilt 

(3.18) u* := uy (y*) + Kev* eB, 


und u*, y* geniigt den Operatorgleichungen (3.16). 


Beweis. Die Behauptung folgt aus (3.17) unter Beachtung von (IV’), (I’), (V’) 
und der aus (II’) und (III’) folgenden Bezeichnung K = K} K}. 


Bemerkung 2. Unter der Voraussetzung 1 (Voraussetzung 2) gelte fiir die L6- 
sung “*, y* von (3.16) 


~ 


(3.19) M,(u*)EB,, KM,(u*)eD(L,). 
Dann ist u* eine ,,klassische‘‘ Lésung von (3.1). 
Beweis. Nach Voraussetzung gilt 
u* = uo(y*) +K M, (u*) 


oe Ly ug(y*) = M,u*. 


Nun ist nach (IV) bzw. (IV’) u)(y*)€D(L,) und es gilt 
(L,u%(y*))(s)=0 fiir seQ. 
Aus (3.20) folgt daher mittels (3.19) unter Beachtung von (II) bzw. (II’) u* €D(L,) 
mit 
(L, u*) (s) = (LZ, K M, (u*) )(s) =M,(w*)(s) (se Q), 
Lg u* = Ly (y*) +LyK M, (u*) = Lot (y*) = My ut. 


Damit ist Bemerkung 2 bewiesen. 


Definition 1. Sei u*, y* mit u*€Y? bzw. u*€B,, y*eB Losung von (3.16), 
dann heiBt w* eine ,,verallgemeinerte‘‘ Lésung von (3.1). 


Aus dieser Definition und Bemerkung 2 folgt 


Lemma 3. Unter der Voraussetzung 1 bzw. Voraussetzung 2 sei u* eine ,,ver- 
allgemeinerte‘‘ Losung von (3.1) und es gelte (3.19). Dann ist w* eine ,,klassische 
Lésung von (3.1). 


Im folgenden werden daher nur Existenzsatze fiir ,,verallgemeinerte‘’ Lésun- 
gen bewiesen. Es gilt der 

Satz 1. Die Voraussetzung 1 sei erfiillt. Dann gibt es mindestens eine ,,ver- 
allgemeinerte‘‘ Losung von (3.1). 


Beweis. Nach Lemma geniigt es die Existenz einer Lésung v*€¥?, y*eB 
von (3.14) nachzuweisen. Wir wenden Hilfssatz 1 an und setzen B= Bb, H = Y?, 


R=b|L|p.5, S(v, y) =v —H#* M, (u(y) +30), T(v, y) = LM, (u(y) +2). 
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Uberpriifung der Voraussetzungen des Hilfssatzes: 
(Aa) T(-,-) ist nach (III), (IV), (I) und (V) stetig von #? x B in B. 
(AB, y) S(-,-) ist nach (III), (IV) und (I) stetig von #? x B in #?, woraus 
(Af) und (Ay) folgen. 
(B) Sei y€B und seien v,, »,€¢ Y?, dann folgt aus (VI) 
(S (v1, ¥) —S (v2, Y), % —v,) = ||, — vq 
—(#* pie (Uo (y) +H 0,) — M, (u(y) +H0,)), Yy — v,) = |v, — rele 
— JX G4 y) +24) — M, (ig (y) +222) (8), (#(v%, —22))(s)> as 
= I = an —a | (v, — 22) | = (1 — 4 | [. e) |. — Vel 
2C |r, — rel, 
womit Voraussetzung (B) bewiesen ist. 
(C) Sei y€B mit ||y|,<R, dann folgt aus 
SO~a;, y) =H* M, (uo (¥)) 
nach (IV), (1) und (III) Voraussetzung (C). 
(D) Sei ve #? und y€B, dann folgt nach (VII) 
IT@, ») le S [Loz [Me (uo (y) +2) |p Sd |Lp.2= R. 


(E) Sei ve F?, a mit |lv||4.< A, |ylzg <A, dann existiert eine Konstante 
A,>0, so daB [u(y )+H rl gS <A, gilt. Nach Voraussetzung (VIII) ist dann 
T(-, -) kompakt von #? x Bin B. 

Die Voraussetzungen von Hilfssatz 1 sind daher erfiillt. Es existiert daher 
mindestens ein v*€Y?, y*€B von (3.14), womit Satz 1 bewiesen ist. 


Ein analoger Satz gilt unter der Voraussetzung 2. 


Satz 2. Es sei Voraussetzung 2 erfiillt. Dann gibt es mindestens eine ,,ver- 
allgemeinerte‘‘ Losung von (3.1). 


Beweis. Nach Lemma 2 geniigt es die Existenz einer Lésung v* € #2, y*€B von 
(3.17) nachzuweisen. Dazu wendet man wieder Hilfssatz 1 an und setzt B= B, 
H=2Z*, R 0 [Llosa S(v, y) =v —K}M, (u(y) + Kte), T(v, y) =LM, (049 (y) + 
Kjv). Die Uberpriifung der Voraussetzungen des Hilfssatzes 1 verlauft ganz 
analog wie beim Beweis von Satz 1. 


Um Existenz- und Eindeutigkeitssatze fiir ,,verallgemeinerte‘’ Lésungen von 
(3.1) aufzustellen, bendtigt man weitere Voraussetzungen. 


Voraussetzung 3. AuBer Voraussetzung 1 gelte: 


(IX) Es existieren Konstanten Rj >0 und Mz >0, so daB fiir yeB mit 
vz S Ro gelte 


4, (1g (¥) )Ilva <= MR||¥*|_e0.g- 


~ 
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(X) Es existiere eine Konstante L Rr, > 0, So daB fiir ve, y,EB (j =1, 2) mit 
lee SM, /C, |y;1p S Ro gelte 
|, (#0 + U9 (¥4)) —M, (Hv + U9 (V9) )Ikga 
= Lp, ||%o (1) — Mp (V2) llve/ (| * oases I“ols+r) : 


(XI) Es existiere eine Konstante k > 0, so daB fiir v,€L%, y,€B (j =1, 2) mit 
lrile SMazJ/C, |lvjlp < Ro gelte 


| (4, + Ug (14) —M, (#05 + Up (V2) ) lle 
s RI(|F loge Lee ae) | (vy — Vg) + Up (4 — V2) lop 


mit 
k (Lx,/ C + [Ito lla! |# lege) <1. 
(XII) Aus Kz=Og» mit z€ YL" folge z=O gu. 
Voraussetzung 4. AuBer Voraussetzung 2 gelte: 


(IX’) Es existieren Konstanten R,>0O und Mz, >0, so daB fiir yeB mit 
v\lz S Ry gelte 


94, (“0 (¥))lle, SMe! |Kb ls, 


(X’) Es existiere eine Konstante Lp >0, so daB fiir ve Y?, y,€B (j =1, 2) 
mit |v: MR /C, ||y;le S Ry gelte 


|, (Kivu + U9 (¥4)) —M, (Kyu + Ug (v2) )Ila, 
= Le, | (V1) — Mo (V2) le,/ (|AGla,.¢ “ole+s,) : 


(XI’) Es existiere eine Konstante k > 0, so daB fiir v,eY?, y,€B (j =1, 2) mit 
lrjleo SMR IC, |lyjle SR gelte 


|, (Ki V1 + Up (14) —M, (Ki Up + Up (Yo) Yile 
S/ (Kila, |Zle2) A$ (v, — v2) + 9 (1 — Ye) le, 


k (Lr J|C +| Uol|5-+2,/||Kile-+z,) Saekes 
(XII’) Aus Kz =O, mit z€B, folgez =O, . 


mit 


Es gilt der folgende lokale Existenz- und Pseudo-Eindeutigkeitssatz: 


Satz 3. Es sei die Voraussetzung 3 erfiillt. Dann gibt es genau eine ,,ver- 
allgemeinerte*’ Losung 


(3.21) Uu* = Uy (y*) +#v¥ € LP 


von (3.1) mit y*€B, v*¢Y?, fiir welche ||y*||, < Ry gilt. Ferner konvergiert das 
folgende Iterationsverfahren 


Vyiy = A*M, (40,44 + U9 (y,)) 


ses AS 
ze) Vy41 = LM, (9,41 + up (¥,)) (v=0, 1,2...) 


mit y,€B und |/yol, < Ry gegen die Elemente v*, y*. 
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Beweis. Wir betrachten zunachst die Gln. (3.14) und beweisen die Existenz von 
genau einer Lésung v*, y* mit ||y*|,<R, und die Konvergenz des Iterations- 
verfahrens (3.22). Dies folgt durch Ravens von Hilfssatz 2. Man setze mit den 
Bezeichnungen des Hilfssatzes: 


S(v, y) =v — #*M,(Hv+u(y)), Tv, y)= LM, (#0 +u9(y)). 
Die Bedingungen (A’«, f), (B’) und (D’) folgen analog wie beim Beweis von Satz 1. 
Uberpriifung der restlichen Bedingungen des Hilfssatzes: 
(A’y) Sei ve F?, y,€B (j =1, 2) mit |o||y.< M/C und |y,|, < Ro, dann folgt 
mittels (X) 
|So(v, V1) —So(% V2) lee = |* (™, (Hv + tg (V1)) — My (Hv + U9 (¥2))) lhe 
= | * loa Le, 
= Lr, ba c= Valle . 
(C’) Sei y€B mit |y|z< Ry, dann gilt nach (IX) 
|S (Og, ») gs = |2P* My (0 (¥) lex SMa. 


(E’) Sei v,¢2?, y,eB (7 =1, 2) mit lojleo SMR IC, 
man mittels (XI) 


(1 —%e )I.ee! ( | "lng \\“olp+oe) 


lz = Ry, dann erhalt 


k 
|Z (,, V4) —T (v2, Vo) le = lace ea Ug) + U9 (V1 — Ve )ILee 
<hr, —%ally + (F “ols+vel|4 ae) ll — Yelle 


? 


womit mit 
h=k, ky=h|\uols,orl|4 lege 
Bedingung (E’) erfiillt ist. 

Es existiert daher nach Hilfssatz 2 genau eine Losung v*, y* von (3.14) mit 
y*|lz Ro, welche nach (3.22) konstruiert werden kann. Nach Lemma 1 folgt 
hieraus die Existenz einer ,,verallgemeinerten“ Lésung w*, die in der Form (3.21) 
mit ||y*||, S Ry darstellbar ist. 

Wir zeigen nun, daB es genau eine enamalleen acute Lésung u* € ¥? von (3.1) 
mit ||y*|, SR, gibt. Dazu sei angenommen, daB u¥*, y* (j =1, 2) mit || y*|,< Ry 
zwei Lésungen von (3.16) seien, d.h. es gilt wegen (V) 

(2 23) uj ar Uy (V;*) +KM, (u;‘) 

' yf =LM,ue. 

Aus der ersten dieser Beziehungen folgt fiir u* die Darstellung 
(2.24) uj mead (v7) + Kz 

mit z*€Y%. Nach (XII) folgt aus (2.23) und (2.24) die Beziehung 


Wendet man hierauf den Operator #* an und setzt 


* __ yx o* 
Uj; =H* 7H, 
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so erhalt man wegen (2.24) und (IIT) 
vf =H* M, (HF + U9 (y;*)) 
y# = LM, (v¥ + u9(y*)). 


v;*, ¥* (7 =4, 2) ist also Lésung von (3.14), und nach dem ersten Teil des Beweises 
existiert genau ein v*, y* von (3.14) mit ||y*||,; S Ro. Es gilt daher v* = v**, y* = y* 
und somit u* = u* (j =4, 2). 

Damit ist Satz 3 bewiesen. 


Ein analoger Satz gilt unter der Voraussetzung 4. 


Satz 4. Es sei die Voraussetzung 4 erfiillt. Dann gibt es genau eine ,,ver- 
allgemeinerte“’ Lésung 


(3.26) u* = Uy (y*) + Ki v*¥ eB, 


von (3.1) mit y*eB, v*¢Y?, fiir welche ||y*|, < Ry gilt. Ferner konvergiert das 
folgende Iterationsverfahren 

Us Ki M, (Ki Uy 44 + Uo (y,)) 

V1 = LM, (Ay Vpti + Mo (y,)) 

mit yy€B und ||/yollz S Rp gegen die Elemente v*, y*. 

Beweis. Man betrachtet zunachst die Gln. (3.17). Der Beweis verlauft dann 
analog wie derjenige von Satz 3. 

Bemerkung 3. Sind die Bedingungen (IX)—(XI) von Voraussetzung 3 bzw. 
([X’)—(XI’) von Voraussetzung 4 fiir alle Ry = R* mit einem geeigneten R* >0 
erfiillbar, dann besitzt das Problem (3.1) genau eine ,,verallgemeinerte“ Lésung. 

Damit diese Voraussetzung erfiillt sein kann, muB notwendigerweise eine der 
beiden folgenden Bedingungen gelten: 

(a) Es ist R=O, dh. es existiert /€B, so daB fiir alle we Y? bzw. we B, gilt 
M, u=f (lineare Randbedingung). 

(b) Es existiert eine Konstante L* > 0, so daB fiir alle Ry = R* gilt Lp SL*, 
d.h. M, geniigt einer globalen Lipschitzbedingung. 

Beweis. Die Behauptung folgt unmittelbar aus Satz 3 bzw. 4 unter Beachtung 
von Voraussetzung (XI) bzw. (XI’). 


4. Anwendungen 
In diesem Abschnitt sollen zwei Anwendungen der allgemeinen Existenzsatze 
gebracht werden. 
Beispiel 1. Als erste Anwendung betrachten wir das nichtlineare Randwert- 
problem 
—Au(s)=f(s, u(s)) fiir s€Q 


ae u(s)=(M,u)(s) fir seaQ. 
Hierbei bezeichne Q eine offene beschrankte konvexe Menge aus R?, 02 den Rand 
von 2 und A den zweidimensionalen Laplace-Operator. 


33* 
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Wir machen nun die folgende Voraussetzung : 

(a) Q sei eine offene beschrankte konvexe Menge aus R?, und der Rand 02 sei 
doppelpunktfrei und stetig gekriimmt°. 

Bemerkung 4 (vgl. z.B. [9], S.176ff.). Unter der Voraussetzung (a) gilt fiir 
EC [AQ] mit der Abkiirzung 


(4.2) W,(s) = af if cos 2 if m4) w(t) dt ® 
an 


die Relation 

(4.3) AW,(s)=0 fir seQ, 
und fiir s,€Q, s€dQ und s,—s gilt 

(4.4) = W,, (S,) = #(s) + W,, (3) - 


Ferner ist die Funktion 
1 cos(|s, ¢|, 22) 


A) 


stetig auf 02 x 02 und das Problem 
H(s) +f A (6.4) wl dt=0 


mit w~E€C° [dQ] besitzt nur die triviale Lésung u(t) =0 auf 02. 


Bemerkung 5 (vgl. z.B. [2] Bd. I, S.316 und Bd. II, S. 282ff.). Unter der 
Voraussetzung (a) gibt es eine Funktion (Greensche Funktion) 


(4.5) K(s,t) = log eal + Ky(s, 2) 


mit K,€C°(Q x Q), A, Ko(s, t) =0 fiir s,t€Q, K(s, t)=0 fiir s€dQ, t€Q und 
K(s, t) = K(é, s) fir s +¢ und s, t€Q. Fiir f€C!(Q] folgen aus 


(4.6) u(s)= f K(s, t) f()dt 
die Beziehungen 
—Au(s)=f(s) fiir s€Q 


4. 
ey?) u(s) =0 fir seoQ.? 


Ferner ist die Greensche Funktion K(s, #) positiv definit (vgl. z.B. [2] Bd. I, 
S. 320). 


Es seien ferner die beiden folgenden Voraussetzungen erfiillt: 
(b) f(s, w) sei eine in H = {(s, w)|s€Q, we R}} stetige reelle Funktion mit 


(4.8) (f (s, m4) —f(s, Up) ) (u, — Uy) Sa(u, —Uy)* 


5 Diese Bedingung kann abgeschwacht werden. 

6 Is, t| bedeutet hierbei den euklidischen Abstand zwischen s und ¢ in R? und Ny 
die innere Normale an 02 im Punkte tf. 

7 u(s)=0 fiir se 02 bedeutet hierbei w(s,) +0 fiir s,€Q und s,>se éQ. 
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fiir s€Q, w,¢ R* (j =1, 2) und mit 
(4.9) aSo  Cr=4 alt= 0. 
Hierbei sei A, der kleinste Eigenwert von 

P= AL KG.) 9 at 
mit pe€C*[Q]. 


Ferner existieren Konstanten D,; = 0 (j =1, 2) und 2S < oo, so daB fiir se Q, 
u€ R? gelte 


(4.10) |f(s, v)| SD, +D,|u|?7?. 


(c) M, sei eine (nichtlineare) vollstetige Abbildung von L? [Q] in C° [dQ]. Es 
existiere eine Konstante b> 0, so daB fiir alle wEL? [Q] gelte 


(4.11) || 2¢\|¢0 a0) Sd- 


Wir beweisen nun als Anwendung von Satz 1 das folgende 


Lemma 4. Es seien die Voraussetzungen (a), (b) und (c) erfiillt. Dann besitzt 
das Problem (4.1) mindestens eine ,,verallgemeinerte‘‘ Lésung. 


Bewets. Wir wenden Satz 1 an. Mit den Bezeichnungen des Satzes setzen wir: 


(442a) B:=C[@Q), PP:=L?[(Q], $4:=L1[Q] — 
= (L,u)(s):=—Au(s), (M,u)(s):=f(s, u(s)) (s€Q), 
a (L,%) (s):= lim u(s,) fir s,€2, s€aQ 
mit 
D (Ly) = C?[Q]0C°(Q)], D(L,) := C°[Q], 
pe) B:= 00}, preeual 
(4.12d) Up (¥) (s) := W,(s) (s€Q). 
(4.12e) (K u) (s) pet K(s, thu(t)dt (s€Q). 


Mit diesen Festlegungen iiberpriifen wir die Voraussetzungen von Satz 1. 


(I) Nach (4.12a, b,c) ist wegen |Q|<oo L, ein linearer Operator mit 
DiEjycFH” und. R(L,)c C2) =B. Ebenso ist L, ein linearer Operator mit 
D(L,)¢<D(L,)¢cY? und R(L,)<C® [@Q] = B. M, ist nach (12b) und Voraussetzung 
(4.10) eine stetige beschrankte Abbildung von L?(Q] in L*(Q] (vgl. z.B. [10], 
S. 154 Satz 19.1 Nr. 2) und M, ist nach Voraussetzung (c) eine stetige beschrankte 
Abbildung von L? [Q] in C®[@Q). 

(II) Der erste und letzte Teil der Voraussetzung folgen nach (4.12e) aus 
Bemerkung 5. Es ist noch zu zeigen, daB K ein linearer beschrankter Operator 
von L?[Q] in L? [Q] ist. Dies wird zusammen mit Voraussetzung (III) bewiesen. 
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(III) Es gentigt wegen |2| < co mit einem ¢ > 0 zu zeigen 


[S| Ks, )|Pt?dtds <x 
22 


(vgl. z.B. [7], S. 56 und 52), was unmittelbar aus (4.5) folgt. Hiermit ist Voraus- 
setzung (III) und der Rest von Voraussetzung (II) bewiesen. 


(IV) Bedingung (3.4) folgt aus (4.3), (4.4) und (4.12d). Fir yeC?[@Q] gilt 
nach (4.12) 
|| 49 (¥) op = (J | Lal, t) y(t) at? ds) Ip 


<(f(,f 14.4] aefPas)” |p 


C°(a2)- 


Da nach Voraussetzung (a) die Menge 2 konvex ist8, gilt fiir s€Q, t¢€dQ 


1 cos(|s, ¢|, 7%) 
It |s, ¢| 


AG >0 


? 


und wegen (4.4) und 
[A(s:2) dia oftirgssecQ 
an 

(vgl. z.B. [6], S. 130) erhalt man 


LAG dD) dta wehiir. scQ- 


02 


Zusammenfassend ergibt dies die Abschatzung 
Move S(f (JA, 9 at's)" Ly 


<2| 2)" Lycra) 


C°[8Q] 


womit Bedingung (IV) bewiesen ist. 
(V) Zunachst folgt mit y€C° [02] aus (4.4) und (4.12b, d) fiir s€¢aQ 


(L2%o(¥))(s) = y(s) + J Als, t) y(t) dt=:(I +A) y. 


Nun ist nach Bemerkung 4 die Abbildung A vollstetig auf C°(@éQ] und —1 ist 
kein Eigenwert. Nach dem Alternativsatz fiir vollstetige Operatoren existiert 
daher eine beschrankte Inverse L : = (I + A) auf C° [dQ], d.h. 


ho ae a | (I +-A)*}p 53 < oO, 


womit Bedingung (V) bewiesen ist. 
(VI) Sei u,eL?[Q] (j=4, 2), dann folgt wegen p=2 und f(s, u,(s))€L7(Q] 
(vgl. (I)) aus (4.8) 


a & (M, (4) —M, (145) ) (s), (%4 — Ug) (s)> ds 
aa (f(s, u,(s)) —f(s, Ue (s))) (u (s) —Ug (s))ds Sa 20, — ell. 


8 Die Voraussetzung, dafS 2 konvex ist kann abgeschwacht werden. 
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Die zweite Bedingung von (3.7) folgt aus (4.9) unter Beachtung von 


1 
||. 19) = —=- 
| I [(Q]L?(Q] VA, 

(VII) Die Bedingung folgt aus (4.11). 
(VIII) Da nach Voraussetzung (c) M, vollstetig ist von L? [Q] in C° [dQ] und L 
beschrankt auf C°[0Q] (vgl. (V)), ist 2M, kompakt von L? [Q] in C°[@Q]. Damit 


sind die Voraussetzungen von Satz 4 erfiillt. Die Behauptung von Lemma 4 folgt 
aus Satz 1. 


Nach Lemma 4 existieren somit Funktionen y*€C°{02] und u* mit 
J | w*(s)|? ds < cv, so daB gilt 
2 


(4.13) u*(s) = f A(s, t) y*(t)dt+ f K(s, t) f(t, u*()) dt 
aa Q 
fast iiberall auf Q und 
(4.13 b) y* (s) gt A(s, t) y* (t) dt = My (u*) (s) 
a 
identisch fiir s€eQ. 


Um die Existenz einer Lésung von (4.1) im ,,klassischen“ Sinn zu sichern 
machen wir die folgende zusatzliche ees eee 
(d) Es seien i {(S.4) (6 ==45 2), und = f(s, uw) stetige Funktionen in §. 


Fir wel? [Q) gelte M, wEC* [dQ] mit «>0 (Klasse der «-hélderstetigen Funk- 
tionen). 


Dann gilt der 


Satz 5. Es seien die Voraussetzungen (a), (b), (c) und (d) erfiillt. Dann gibt es 
in der Aquivalenzklasse der Lésungsfunktionen aus L? [Q] von (4.13) eine Funktion 
u* €C2([Q\>C?[Q)} mit 6>0, welche das Problem (4.1) im ,,klassischen“ Sinn lost. 


Beweis. Fiir s,, s.€Q folgt wegen f(t, u* (t))€L?[Q] 
A K(s,, é) f(t, w*(t)) dt ih K (sq, t) f(t, u* (t)) ad 
<(f| Ks.) —Klse, Ohad) ( f|t(t w* (t))|? de)" +0 
2 2 


fiir s,;—>s,, d-h. es gilt 


(4.14) J K(s, t) f(t, w* (é)) dt€C°[Q)]. 

2 
Nach (4.13) kann daher wegen (4.3), (4.4), (4.12) und (4.14) aus der Aquivalenz- 
klasse der Lésungsfunktionen von (4.13) w* €C°[Q] gewahlt werden, und mit die- 
sem u* ist die Beziehung (4.13) identisch auf  erfiillt. Aus (4.13) folgt daher 
mittels (4.4), (4.14) und Bemerkung 5 fiir s€Q und s—s¢e0Q 


lim u* (s) = Sam J At (s, t) y* (0) dt, 


Sos 
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d.h. zusammen mit (4.13 b) erhalt man 
u* (s) = M, (u*) (s) 


fiir s€dQ, womit die zweite Beziehung von (4.1) bewiesen ist. Wegen u* €C® (Q] 
ist f(s, u*(s))€C®°[Q], d.h. es gilt (vgl. z.B. [2], Bd. I, S. 317) 


(4.45) if K(s, t) f (¢, w* (t)) dteC?[Q]. 
Ebenso ist wegen y* €C® [dQ] (vgl. z.B. [5], S. 47) 
Zar t) y* (t) dteC* [0Q] 
mit ¢>0 und somit folgt wegen M, (u*) €C* [dQ] (Voraussetzung (d)) aus (4.13 b) 
y* €C® [GQ] mit 6>0. Hieraus ergibt sich (vgl. [5], S. 48) 
Als t) y* (t) dteC*[Q), 
und nach (4.13a) und (4.15) gilt daher u*€C°[(Q]. Nach Voraussetzung (d) ist 


somit f(s, w*(s))€C®[Q] und nach dem Satz von Poisson erhalt man schlieBlich 
(vgl. [5], S. 87) [ K(s, #) f(t, u* (t))dt€C?[Q] mit 
9 


—A f K(s, t) f(t, w*(t))dt=f(s, w*(s)) fir seQ. 
2 
Aus (4.13a) folgt daher mittels (4.3) die erste Beziehung von (4.1). Damit ist 
Satz 5 bewiesen. 
Bemerkung 6a. Es seien k(s, 2), ae k(s, t)€C°[8Q x Q] und a(s), b(s) €C1[2Q]. 
Fir w€L? [Q] wird beispielsweise durch 


(M,w) (s) := a(s) +(s) /(1 ip (J h(s, t) u(t) dt)’) (s €€Q) 


eine Abbildung M, erzeugt, die den Voraussetzungen (c) und (d) geniigt. 
Bemerkung 6b. Es sei c(s), d(s) €C![Q] mit d(s) =O auf Q. Ferner sei 


f(s, u) =e(s) —d(s) w'"—*/(1 +0) 


mit m=2 und ganzzahlig. Mit p=2m, a=0 und C=1 sind dann die Voraus- 
setzungen (b) und (d) erfiillt. 

Beweis. Die Bedingung (4.8) folgt mittels des Mittelwertsatzes der Differential- 
rechnung. Der Rest ist trivial. 


Als nachste Anwendung betrachten wir das 
Beispiel 2. Vorgelegt sei die nichtlineare Randwertaufgabe 


—u''(s)=g(s,u(s)) fir seQ:=)o, 1[, 


se u()=K@), u(t) = Bw). 
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Wir machen nun die folgenden Voraussetzungen: 
(a’) g(s, #) sei eine in Hy) = {(s, u)|s€Q, wER1 stetige reelle Funktion mit 


(4.17) (g(s, U) —g(s, Uy) ) (Uy — Uy) Sa (uy — Up)? 

fiir s€Q, u;€R} (j =4, 2) und mit 

(4.18) a=0, Ci=1—— > >0. 

Ferner existieren auf Q stetige Funktionen g,(s) (j=14, 2) mit f (g,(s))?ds< co 
Q 


und eine auf R} ={w|weR!, w=0} stetige Funktion h(w), so daB fiir s€Q und 
u€R? gelte 


(4.19) |g (s, #)| Sgr (s) +82(s) (|). 

(b’) Es seien F; (j =4, 2) stetige Abbildungen von C°(Q] in R} mit 
(4.20) |B (u)| Sb (7=4, 2) 
fiir alle w€C°[Q] mit einer geeigneten Konstanten b> 0. 


Als Anwendung von Satz 2 beweisen wir das 


Lemma 5. Die Voraussetzungen (a’) und (b’) seien erfiillt. Dann besitzt das 
Problem (4.16) mindestens eine ,,verallgemeinerte‘‘ Losung. 


Bewets. Wir wenden Satz 2 an. Mit den Bezeichnungen des Satzes setzen wir: 


BP Be= TO} Be= C10 “BER 


ee) mt dex Norm OPE tar’ gee 
ney. (L, ¥) (s) = —u'"(s), (M, ) (s) =g(s, #(s)) (s€2), 
ete) L,u = {u(0), w(4)}, M,u={Fu, Ru} 
mit e 
by D(L,) =C?[Q]0€°[Q], D(L,) =C°[Q], 
eae B= CQ), By = C°(Q]0L2(0}, 
(4.21 d) Uy (y) (Ss) := Vy + Ves fir séQ und y= {y,, ya}eR?, 
(4.21 e) (K wu) (s) ee thu(t)dt fiir seQ 
mit 
(4—s)¢#, OStSs 
(4.21) aie en 


Mit diesen Festlegungen iiberpriifen wir die Voraussetzungen von Satz 2. 

(I’) L, ist nach (4.21, b, c) ein linearer Operator mit D(L,)¢ B,¢L?[Q] und 
R(L,)¢C°[Q]=B. Ebenso ist L, nach (4.21a, b,c) ein linearer Operator mit 
D(L,)¢D(L,) = B, und R(L,)¢B. Da nach Voraussetzung (a’) g(s, u) stetig ist 
in §) und Bedingung (4.19) gilt, ist M, nach dem Satz von Lebesgue stetig und 
beschrankt von C°(Q] in L?[Q]. Nach Voraussetzung (b’) ist M, ein stetiger 
beschrankter Operator von C°[Q] in B. 
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(II’) Sei f¢ By, und fiir s€Q gelte 
=f K(s,t) fd dt, 
2 
dann folgt durch elementare Rechnung 
—u''(s)=f(s) (s€Q), u(0)=u(1)=0. 
Sei v€L?(Q], dann gilt fiir s€Q die Darstellung 


(4.22) (Kv) (s) := >) Pe) o,(s) 
i=1 
mit 
(4.22b) pi(s)=y2sn ais, A= Pa. 


Der durch (4.21 e, f) erzeugte eens K ist linear und beschrankt von L?/Q] 
in C°[Q]. Ebenso wird nach (4.22) durch K (s, ¢) ein linearer, beschrankter, selbst- 
adjungierter, positiv definiter Operator K auf L?(Q] erzeugt, fiir welchen die 
Beziehung J, K =K besteht, (J, ist die Injektionsabbildung von C® [Q] in L?[Q)). 
(III’) Da J, die identische Abbildung ist, gilt K} = K+ (eindeutig bestimmter 
positiv definiter Wurzeloperator von K). Fiir v€L?[Q] setzen wir mit seQ 


edu v) 


t=1 


Es bleibt noch zu zeigen, da8 A} linear und beschrankt von L2[Q] in C°[Q] ist. 
Mittels der Schwarzschen Ungleichung erhalt man 


/ oo 4 co 
M h < " 2) M Gi ( (s-Eh) — gi (s s))? 
us Bia dea Us! veils (De 2 ) ates, ee ms 
2)2 foe) sin?i > h\8 
= = 2 72 7 ele: [2]° 


Die Reihe > o; (h) mit s,; (A) = “s sin? 7 = h ist absolut und gleichmaBig konver- 
i=1 
gent, und fiir jeden einzelnen Summanden gilt s;(h) +0 ({=1, 2,3...) fiir h->0. 


Durch Vertauschung der Reihenfolge der Grenzwerte folgt daher 


h—0 s,st+heQ 14 


t 
lim Max_|(K}v)(s +4) —(K}v)(s)| < 9322? (3 tin lim s ) loll ,o] = 0- 


Es gilt daher Kj v€C°(Q]=B,. Mittels der Schwarzschen Ungleichung erhilt 


man ferner 
$ co $ 
(pi (s))? 
D> Max Ae 
of s(S (9% ») seQ (> Ai 


PSS V|\¢° (oles Sn 


< (pi, ¥ pe 
pap ie 


V2 1 
2 (5 seu, Pleo 


141= 


womit (III’) vollstandig bewiesen ist. 
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(IV’) trivial. 


(V’) Eine elementare Rechnung zeigt, daB L = fa i} die Bedingungen von 
Voraussetzung (V’) erfiillt. 


(VI’) Sei u4,€C°(Q] (j =1, 2), dann folgt wegen (4.17) 
CAVACA) — JM, (ue), Ait — Siz) = of (g(s, #1 (s)) —g(s, ty (s))) 


. (m4 (s) —w(s))ds< a (4 (s) — te (s))?ds Sa|J, (u —%2) lis (ay) 


womit (3.11) bewiesen ist. (3.12) folgt aus (4.18) unter Beachtung von 


eet 


= oi Vi, eet 


(VII’) Sei w€C°[Q], dann gilt nach (4.20) 


M, (#) || = Max (|K(v)|,|B(#)|) So. 


(VIII’) LM, ist kompakt von C°[Q] in B = R?, da R? endlich dimensional ist. 
Die Voraussetzungen von Satz 2 sind erfiillt. Es existiert daher mindestens eine 
, verallgemeinerte’’ Losung von (4.16). Damit ist Lemma 5 bewiesen. 


Nach Definition 1 und Lemma 5 existieren daher ein y* = {yf, yf}¢R? und 
ein u* €C°[Q] mit 
ft 
(4.23 a) u*(s)= yi +yes+ f K(s, 2) g(é, w*(f)) de, 
0 
(4.23 b) v=Ku*), vi tyt = Blu"). 


Wegen w*€C°[Q] und Voraussetzung (a’) gilt g(s, w*(s))€C°(Q]>L?(Q] =B, 
Aus (II’) folgt hieraus KM, (u*)€D(L,). Zusammen mit Bemerkung 2 ist u*(s) 
eine ,,klassische‘‘ Lésung von (4.16). 


Zusammenfassend erhalten wir den 


Satz 5. Es seien die Voraussetzungen (a’) und (b’) erfiillt. Dann gibt es min- 
destens eine Funktion w* (s) €C®[Q]>C?[Q], die das Problem (4.16) (im klassi- 
schen Sinn) lost. 


Bemerkung 7a. Es seien a, %, .--, %, reelle Zahlen, dann erfiillen beispiels- 
weise die Abbildungen 


as 


1 , 
1 + ag (w(0))?+-09 exp (w(E)) Hos f (u(s))?ds 


F (vu) =a, + 


Oh (U (0))® tag 


iF al (eV Faa (1)? bags Max |1«(s)|_ 


F, (u) = % + 


die Voraussetzung (b’). 
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Bemerkung 7b. Es seien a(s), b(s) und c(s) auf [0, 1] stetige Funktionen mit 
b(s) =0 und c(s) =O. Ferner seien die Konstanten B;< } (1 =1, 2, 3, 4), und es sei 


b(s) ; c(s) yrzm—-t 
6 (Si) = te ie yeas iar ae A 


mit m=2 und ganzzahlig. Dann ist die Voraussetzung (a’) erfiillt. 


Bemerkung 8. Die Probleme unter Bemerkung 6 (mit d(s)=0 und m= 3) und 
Bemerkung 7 (mit 6(s)=:0 oder c(s)==0 und m=3) kénnen mit Hilfe der Exi- 
stenzsatze von Ehrmann [3], [4] und Conti [1] nicht behandelt werden. 
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Abstract. Approximation of a conformal mapping of a simply connected domain 
onto a circle by means of Ritz’s method yields a system of linear equations with a 
Gram matrix. The asymptotic behaviour of the minimal condition of this matrix is 
studied in dependence on its order. 


i: 


Introduction. The problem to determine numerically a conformal mapping 
of a bounded simply connected domain G onto a circle can, with the aid of 
extremal functions, be reduced to the solution of a sequence of finite systems 
of linear algebraic equations with Gram matrices A, (n=1, 2,...), the order 
of A, being ~. The elements of these matrices depend upon the definition of 
the Hilbert space H where the solution is sought and upon the basis chosen in H. 

The influence of rounding errors on the numerical solution of the system 
A,,x = 6 (especially for Gauss elimination in floating-point computation) may be 
studied in terms of Bauer’s scaled condition number 


(4) B(A,) = inf cond (D, A,, D,) 


where cond(B)=|B| |B and D,, D, are diagonal matrices of order » with 
positive diagonal elements (for real matrices see [1, 4]). By || Bl] we shall denote 
the spectral norm of B. In order to maintain the Hermitian symmetry of the 
matrix, we shall moreover suppose D,=D,. This is in fact no restriction, since 
the value of the infimum in (1) does not change (for proof see [3]). 

The aim of this paper is to contribute to the solution of the problem if there 
exists a reasonably computable basis of a given Hilbert space of analytic functions 
satisfying the inequality 

sup B(A,,) << ~. 


A positive answer to this question is obvious if G is a disk and H=L,(G). How- 
ever, as we shall show, this inequality does not hold for domains G with boundaries 
arbitrarily close to the circle, e.g. ellipses and polygons, with the same definition 
of base functions. 


The Variational Method. Let G be a bounded simply connected domain in 
the complex plane, 0€G. Let w=f(z) be the conformal mapping of G onto a 
disk with center at the coordinate origin, satisfying the conditions /(0) =0, 
(0) =1. (Then the radius R of the disk is uniquely determined.) 
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Let L,(G) denote the set of all the functions F which are regular in G and 
satisfy (using notation z= x +77) 


\FP= ff |Fe|taxdy < oo. 
If F,, F,€L,(G), we can define the scalar product 
(A, Fs) we! SF, (2) F,(2) dxdy. 


By this definition the set L,(G) becomes a Hilbert space (see [5], p. 117). 

It follows from the least area principle of the theory of functions (see [2], 
p. 189, [5], p.119, [6], p. 281) that the derivative /’(z) is the only solution to 
the problem of minimization of the functional 


JS|F@) (2)|?dxdy 


on the set of functions FEL, (G) satisfying F(0) = 


The derivative /’(z) can be approximated by means of the Ritz’s method 
({5], p. 120; [6], p. 383). We shall give here a slightly more general formulation 
of this method than the mentioned monographs do. However, it is easy to see 
that the generalization does not influence the validity of the assertions cited 
without Hie 

Let {uy n(2)} (k=O, 1,...,; m=1, 2, ...) be a complete system of functions 
in L,(G), such that for ana m the functions u, ,,(z) (ke =0, 1, ..., 2) are linearly 
independent. Moreover, let u)(0) ==0. Then the -th approximation of the deriv- 
ative /’(z) is defined to be the function ¢,(z) which minimizes the functional ||F|j 
in the class of all complex linear combinations of the functions uw, u,,..., u 
satisfying F(0) =1. The function q,, is uniquely determined by the relation 


n 


(Pro Yn) =O 


which is satisfied by any linear combination yp, of the functions w, ,, fulfilling 
the condition yw, (0) =0. If we write 


n 
Pn = OS Aen UE n> 
k=0 
then, under additional assumptions 
Uo (0) =1, Ux n(0) =0 eS Ree Ree) 


the coefficients a,,, (k=1,2,...,) are given by the solution of a system of 
n linear equations ; 


n 
(2) 2 (Up, ni Un ) ay, ae = Wee U; n) qG= 1.72, o0s5 n) 


with a Gram matrix A,,= ((#, ,, 4; ,))- Hence A, is a positive definite Hermitian 
matrix for each n. 

In practical eae eli it is common to use the polynomial basis of L(G), 
ie. to take u,(z)=z* (k=0,1,...). (We have dropped the subscript , since 
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the base functions do not depend on n.) If G is a disk with center at the origin 
and radius R, then the matrix A, of system (2), generated by the polynomial 
basis, is a diagonal matrix with diagonal elements R?“*”/(k +4) (k =4, 2, ..., 2) 
hence (denoting the maximal and minimal eigenvalue of A, by A,,.(4 
Amin(A,), Tespectively ) 


, 


) and 


max n 


lim cond (4,,) = lim Amax (An) a oy 
eo n—>oo Amin (A a) 


However, if we prefer the base functions v,(z)=C,z* with C,=V(k+1)/R* 
(k=0, 1,...), then, denoting B,, = ((v,, v,;)), we have 

cond (B,,) =1 
for each n. 


Obviously, the replacement of the system {u,} by {v,} is equivalent to the 
scaling of A, by the diagonal matrix D,,(C,,..., C,), ie. to the replacement of 
A, by Ds re Eee 

Hence if G is a disk with center at the origin and the matrices A,, are generated 
by the polynomial basis, we have 


sup B(A,)< «~. 


In the following we shall show that this inequality does not hold for a large 
class of domains not only for the polynomial basis, but for a more general class 
of bases uy = 1, iC, (k=1, 2,...,n;n=4, 2,...), where C, ,, are arbitrary 
complex numbers. 
if: 
For further need, we recall the following 


Lemma. Let A be a positive definite Hermitian matrix and B its principal 
submatrix. Then we have 
Amin(A) S Amin (B) ; 
A A) 2 Ler (B) : 
Definition. A domain G is said to be starlike with respect to a point ¢ if each 
ray initiating at ¢ intersects the boundary of G exactly at one point. 


max ( 


Theorem. Let G be a starlike domain with respect to the coordinate origin, 
situated inside the unit circle and such that its boundary touches the unit circle 
exactly at a finite number of points. Let the functions w, ,, have the form C foes 
with arbitrary non-zero complex constants C,, (k=0,1,...,”; m=1, 2,...). 
Let A,, be the matrix of the system (2). Then 


sup &(A,,) = oo 


where &(A,,) is the Bauer’s scaled condition number (1) of the matrix 4,,. 


Proof. a) Notice first that the system u, ,, is complete in consequence of the 
starlikeness of the domain G (see [7], p. 428). 
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Suppose that the theorem does not hold, i.e. that 
(3) sup &(A,) < o. 


The value of &(A,,) obviously does not decrease with increasing n, hence 


sup &(A,) = lim &A,). 


Since A,, is a positive definite Hermitian matrix, we have 
cond (A = |, A, *| = s” ye | A.) Avia (A,,) Z 
The same holds for the scaled matrices D,,A,,D,,. 


We can choose matrices D* (n =1, 2, ...) so that 


dmax(D* A, Dt) 
jim 4(A,) = jim >" (ps4, Ds) * 


Denote by d,,, the k-th element of the principal diagonal of the matrix Df. 
Let a (7, k=1,..., 2) denote the elements of the matrix Df A,, Df, Le. 


of") = ad; n dp, Ore: aCe sine z or 
We can write 


(4) al) =F Vin Ven Pi, k 
where 

(z*, 23) 
(5) Bie= TAL [ay 


is independent of m and 


ya di nen |=". 


The assumption (3) implies the existence of two positive numbers m, M such 
that we have 


(6) MS Ayig (D¥ A, D8) S Angy (D8 A, D8) SM 
for all ». The diagonal elements of the matrix D* A, D* have the form 
of) =|Ye, 0? (k=1,...,%). 
Hence, according to the preceding Lemma and to the relation (4), we have 
(7) m=|y, ,/?SM 
LOA ee ese 


Let B denote the matrix of infinite order with elements Bi, (7, R=, 2, --.) 


defined by (5). We shall prove the existence of a sequence {B;, x} of non-diagonal 
elements of the matrix B, such that 


(8) Jim [Bal = 1. 


Since all the diagonal elements of the matrix B are equal to 1, this will enable 
us to construct a sequence {X,} of second-order principal submatrices of B con- 
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verging to a matrix 


satisfying | |= 1. To the sequence {X,}, we can choose a sequence {Y,} of equally 
situated second-order principal submatrices of matrices D¥ A, ,D3,, where the 
sequence {#,' satisfies the conditions 


p, = max(t);'R;) , 


vy? v 


lim 2, == oo. 


v—> CoO 


Then the relations (4), (7), (8) imply the existence of a subsequence of {Y,} con- 
ing t natrix = 
verging to a matrix \3|2 30,6 


where 6,, d2 are limits of certain subsequences of the sequences {y,,,,} and 
{Yx,, p,}» respectively. Obviously, the minimal eigenvalues of submatrices belonging 
to the mentioned subsequence converge to the minimal eigenvalue of the limiting 
matrix, i.e. to zero. This implies (by the Lemma) the relation 


inf /.,;,(D* A,,D*) =0 


which contradicts the assumption (3). 


b) Using polar coordinates to describe the boundary of the domain G, we 
can write y=r(p) (where gy =arg z, r=|2|). Then we have 


(9 Bis= parqel = Vonen 
where 

4(j+1) (k+1) 

G+re+2)? 


b= 


m 2m ‘ 
/ [r (p) ¥+* +2 cos [(k —7) ¢] eames [r(y) + **2sin [(k—J) p] dp 
Qir == 


22 . 2x \3 
{i ir (p) BU +Vdp- fr (gp) A+D dp 
0 


Denote the difference k —7 by L. In part b of the proof, we will consider L as 
a fixed number and we will find the limit 


lim Re B; j+1- 


joo 


We shall use the following notation. 


Tir (g)]Ecos(L 9) - r(g)]* 40 


P(x) aan o 2a ? 
Su@)¥de 
1 tr(@)]-Ecos(L 9) r(p)* ag 
P(x)=* 


fe ("do 
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It is easy to see that 
(10) (Re B;, 541)? = 9; jar PB (2G +41)) P_(2G +L +1)). 


For 7 -> co, we have 
De jet om ROE 


We shall now study the rest of the expression in (10). 
Let m denote the number of common boundary points of the domain G and 


the unit circle and let p, (v=1, 2,..., m) denote the arguments of these points. 
Thus we have 7(y,) =1 (v=1, 2,...,m) and r(p) <1 if 
pe Uy. 


First, o being an arbitrary positive number, we shall construct a system of 
intervals J,(o, x) (v»=1, 2,...,m; x >0) satisfying the following three conditions. 


(i) p€J,(o, x); 
(ii) if d[J] denotes the length of the interval J, then for any x >0 we have 


o m 

mn SD AU, M503 
(111) f rdp= f r*do 

Jv (6, x) Ju (a, #) 


108, PA, 2)... Mana all 4.0. 
Define’g;(y) Y=1),2—., m) Dy 


g(y)=r(y, +). 


For a moment, consider x fixed. Let uw satisfy the following equality: 


min f [6,(9)I"dy = J" [gy 9)" dy = ae, 2). 


The continuous dependence of the integral on its upper bound of integration 
implies the existence of numbers o, <o/(2m) such that for y= we have 


Fle) Fay= a(o, x). 


Moreover, define o,, by o,,=o/2m. Now we can construct the right-hand sections 
of the sought intervals: 


Jy,(9, %) = Py Y +9,> 


(v=1,..., m). The left-hand sections J, (a, x) can be constructed in a similar 
way, and finally we take 


Je (o, x) a TeX, x) UJ,_(o, x) 


(v=1,..., m). It is easily seen that these intervals satisfy the conditions (i) to (iii). 
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In the following, we shall assume that o is sufficiently small to make the 
intervals mutually disjoint. 


Let 


a(o,*)=fr*dp (v=14,...,m). 
is 


We shall use the following notation: 


J=J(o, )= 0 Jule, ), 
K= K(o, x) = <0, 2%) — J(a, x). 


Then we can write 


S Seale) dp [esl rap 


v=1 Jv 
m a(o, x) T fray 
: J 
(11) P(x) = fade 
K 
tS frag 
aM 


where g,(y~) =7(g)**cos Ly. By the sign + we mean that one of the signs + 
and — is valid on both sides of the equality. We shall also write g instead of g, . 


By the first mean-value theorem, there exist numbers 9g, (a, x) (v=1,..., m) 
(generally different for P. and P) such that 


js rdp=gl9,(e, “Nf r*dp=glg,(a, x)]a(a, x). 


Hence the first term of the numerator in (11) is equal to 


im 2819410, 2)1. 


Let 
M = maxr(¢). 


peK 
We have M <1, hence for each » (vy =1, ..., m) there exists a set J,* ¢ J, of positive 


measure so that 
inf 7(y)>M. 


gpel¥ 


Define ]* by 
fran): 
y=1 
and let w(J*) be the Lebesgue measure of the set /*. If we write 


N= int 7(¢), 


gpey* 


Q = max ¢(¢) 


gpeK 
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(Q is independent of x), we are led to the conclusion 


Sg(y) dy 
K Q (22—a/m) M* 
(12) rag STINT (eee 
Similarly we have 
Sr*dp 
(13) ea <C(c) (M/N)* > 0. 
df 


In the following, P(x) will denote either P(x) or P(x). Let us now take o 
as a fixed number. For each », all the numbers , (a, x) belong to a closed interval 
the length of which does not exceed o. Hence we can find a sequence {x,} and 
numbers w, (a) such that 

jim P(%,) = lim sup fies 
and 
lim Py (a, Xx) 3 Oo, (c) 


nN—> CO 


for y=1,2,..., m. Then by (11), (12), (13) we have 


lim sup P(x) = lim P(x,) = lim {1/m > glos(o, *4))} = AI ¥ ga,(0)). 


x%—> CO 


However, P(x) does not depend on oa, so that we have 


lim sup P(x) = tm > ¢ [lim w, (a)|=1/m 3 g(y,). 


%—> Oo o—>0o = 


In a similar way, we can prove 


m 


lim inf P(x) =1/m > g(y,), 


*%—>0o vel 


whence (since 7(y,) =1) it follows that 


lim P(x) =1]m >) &, (p,) =1/m D>, cos Ly,. 
%—> CO vy=1 p=l 
Returning to (10), we obtain 


(14) Jim |ReB;,¢+11 = 1)/m| $) cos Ly, 


Remark. This already proves the theorem under the ‘additional assumption 
that y,/z is rational for y=1, ..., m. For example, if G is the interior of an ellipse 
with y,=0 and y,=2, we have (setting L = 2) 


fa |ReB; j+2| =1- 


Since |8;,| <1 for any j,k, we thus obtain relation (8) which (by Section a) of 
this proof) proves the theorem. 
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c) Now it is left to prove the existence of a sequence {L,,} of natural numbers, 
such that 


(15) lim cos (L,,y,) =1 


n—> CO 


for y=1,...,m. This will enable us, by combining (14) and (15), to find se 
quences {j,,} and {k,,} satisfying 


im. | Re Bin | = 1 
and hence also 


im. | Binen| ate 1 i 


Let us first assume m= 1 and write simply p instead of y,. We shall show 
that there exists a sequence {L,, K,} of pairs of natural numbers, such that 


(16) lim |L,y —2K,,2| =0. 
Writing w =22/y, we have 
[Ly —2Ka|=|y|-|L—Ko| 


for any L, K. We can construct the sought sequence by the following iterative 
process. 


Let L,, Kg be natural numbers such that 
Ko Laan. 


where 0< |6)|<1. Let us assume 6)>0. For a negative d) we can proceed 
quite analogously. Then for any 7 we have 


7 Kym =J Ly +7 bo. 
Let us define 7) by 79 =[6, ‘]. Hence we have 
fo 1, (Jo +1) do>1. 


If j,6)=1, then we can define K,= 7) Ky, L,= Jol) +1 (n=, 2, ...) to satisfy 
(16). Therefore we shall now assume 7909 <1. We have 
Jo Kyo =Jolo +1 —4,, 


(79 +1) Ko@ = (Jo +1) Lo +144, 
where 
dy =1—190) = 0, = do= (fo +1) 0, —1>0. 


Obviously d,+d,= 6), whence min (d,, d,) = 69/2. Now if d,<d,, we set 
Ky =o Ko, 


L,=Jolo +1, 
Oi —d,. 
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In the other case, i.e. if d, =d,, we set 


Ky = (jo +1) Ko, 
L,= (jo +1) Lo +1, 
Op Ry. 
Thus we come to the relation 


K,o=L, 4.6; | d,)'s 6,/2. 


If we started with a negative 6), we should obtain the same relation. 


In the following step, K,, L, and 6, take the place of Ky, L, and 69, respec- 
tively. By this process, we obtain a sequence {K,,, L,,} such that 


lim |L,y —K,7| S lim (|p| - 69/2") =0 


n—> CO 


which is equivalent to (15) for »=1. 
Now let m>1. Let us assume that there exists a sequence {L,} satisfying 


(17) lim cos L, y, = 1 
for y=1,...,m—1. We claim that (17) is then true for y=1,..., m. 


We shall use the following notation. ¢=y—2ka, k being an integer such 
that ¢€(—2z, 2), and ¢=q(mod 227). 

Let us consider the following two statements. 

(i) There exists a sequence {L,} of natural numbers, such that 


(18) lim max ; Lvl =I 


N— OO v=1,...,m— 


(ii) To each sequence {L,,} enjoying the property (18), there exists a number 
a> 0, such that 


(19) Vel = 4 
for all n. 


If we prove that these two statements contradict each other, it will imply 
that if condition (i) is true, then it holds also if we replace (18) by 


eee eee 
It is easy to see that this is equivalent to (17) for y=14, ad mM. 

To this end, considering the conditions (i) and (ii) as being true, let {L,} 
be a sequence satisfying (18) and let a>0 be the number corresponding to this 
sequence by condition (ii). 

As a consequence of (i), to any natural number k there exists a number N(f), 
so that we have 


—~Y 
at ed | Ln¥,| = ajke 
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for n = N(k). We can assume that 


lim N(k) = ©. 


n-> Co 
Then for u=1, 2,..., 2% we have 
a Ra 
max | uLy cry Yy| Salk. 


v=1,...,.m—1 


Let us construct a sequence 
{,} = Lyay, Ly), 2Ly ys Lye), 2L nis), 3L ays «+> 
where for n=k(k —1)/2+p, w=41, 2,..., 2 we have 
ty = UL): 


Hence to any natural number k there exists a number N(k) such that 


max | \t,y,| Salk 


v=1,..., 


for n> N*(k). On the other side, by (ii), there exists a positive number b <a 
so that 

[tn Wm| > 0 
for all . 


Let c be a natural number such that a/c <b. We shall select those elements 
of the sequence N(k) where k =qc with a natural g and form the following sub- 
sequence of {t,,}: 


{04} = Leys Leecys 2Lnacy» Luaey» 2L acy» 3 Lwiacys +++ 


ie. for n=q(q —1)/2 +m, w=1, 2,...,9, we have v,=yLyy,). 
Obviously, we have 


max _|0,y,| <a/k=al(qe) <b/q 


v=1,....m—1 


starting with a certain m, for any natural g. At the same time, we have 


(20) [On Pm| > 2 
for all 1. 


If the inequalitv (20) holds, then, having in mind the structure of the se- 
quence {v,,}, we realize that the following assertion must be valid: 


(iii) To any natural number q there exists an integer w, such that the in- 
equality Bus 
[4 Ym| > 2 
poids for wa=41, 2, og. 
We will show that no y,, enjoys the property (tii). 
Let @ be an arbitrary number belonging to the interval (b, 27 —b). Let x, N 
be defined by 


(21) x= ([2n/b] +1, N= [(loga —log b)/log 2], 
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square brackets denoting the whole part of the number. Let us assume that 


22) lig|>6 


for 7 =1, 2, ...,% where qx”. 
Let k* = [2a/p]. Hence we have kj <x —1. By assumption (22) we have 


ee pe 
ki<2a—b, (kf +1)p>0. 
Further let 6) = min (d,, d,) where 
d,=2n—kfp, d.=(ko +1)9 
We have d,+d,=9, whence 


Oo = y/2 < (2% — D)/2. 
Let us define ky by 
is if ad<d, 
aie bai) 16 Gey 
Obviously ky < x. 
For simplicity, we shall now assume d, =d,. In the other case, we should 


proceed in an analogous way (only with the difference that points with arguments 
1Rkyp for growing 7 move along the unit circle in the opposite sense). 


Let kf = [22/6], ie. we have 


ur Pe 
h¥ 09> (RE +1) do. 
At the same time 
kf <“—1. 


Let 6,=min(d,, d3) where 


: pe tel tee 
dy= 20% —Rf 09, dy= (RE +1) do- 


Obviously 
6, S65/2. 
Let us define k, by 
& Be ih eae ay 
Ri +1 if dj2dy, 


so that we have k,< x. 


If we now replace 09, ky by 6,, k,, respectively, and repeat the process, we 
can construct sequences (6,1, {k,} such that 


OS One le. 
(23) n ol ze/ 
<= 2% 
for n=1, 2,.... However, the construction of 6, shows that 


J, =|K,,9| 
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where K,, =k, k....k,,. Hence for all m we have 
(24) fOr et ak 


Combining (23) and (21), we see that 6,<6 holds for »=WN. This implies the 
inequality 


— 
| Ky | Sb. 


By (24), we have Ky <x’, where x’ is independent of the choice of g. This 
proves that no number y,, can enjoy the property (iii), and hence the statements 
(i) and (ii) are contradictory. 


Consequently, there exists a sequence of natural numbers {L,,} such that 


lim cos(L,,) =4 


n—> Co 
Tor v1, 2... , M. 
This completes the proof of the theorem. 
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Hinweise fiir die Autoren 


Die Autoren werden freundlichst gebeten, bei der Abfassung ihrer Manuskripte folgende Punkte 
zu beachten: 


Dem Text ist ein englisches Summary voranzustellen, das fiir Kleindruck zu kennzeichnen ist. 


Die verwendeten Symbole miissen so klar bezeichnet sein, daB auch beim Laien kein Zweifel tiber 
Stellung und Deutung auftreten kann. Alle Formeln sind méglichst mit dey Maschine zu schreiben, wobei 
darauf zu achten ist, da8 Indizes und Exponenten trotz des fehlenden GréBenunterschiedes genau als 
solche zu erkennen sind. Anderenfalls miissen sie in geringerer Gré®e mit der Hand eingetragen 
werden. Besondere Lettern (griechische, gotische, Script) sind durch farbige Unterstreichung zu kenn- 
zeichnen. Zur Vermeidung von Verwechslungen wird empfohlen, griechische Buchstaben rot, gotische 
Buchstaben blau und Scriptbuchstaben griin zu unterstreichen. Es wird gebeten, auch gotische und 
Scriptbuchstaben mit der Maschine alslateinische Buchstaben zu schreiben und allein durch die farbige 
Unterstreichung zu kennzeichnen. Sofern dennoch handgeschriebene Buchstaben vorkommen, unter- 
scheide man (auch bei lateinischen Buchstaben) groBe und kleine Buchstaben; groBe Buchstaben sollten 
zweimal, kleine einmal unterstrichen werden. Dies ist besonders wichtig beic, C;k, K; 0,0; p, P; s, S; 
u, U;v, V;w, W; x, X; z, Z. Besonders sorgfaltig beachte man die Schreibweise bei handschriftlichen 
Buchstaben, wenn gleichzeitig e und / oder uv und » oder » und yr oder o und 0 als mathematische 
Bezeichnungen auftreten. Auch gleichzeitig auftretende v und v sowie € und e¢ geben zu Verwechs- 
lungen AnlaB. 


Samtliche Buchstaben in Formeln, Einzelbuchstaben im Text sowie unterstrichene Textstellen 
werden automatisch kursiv gesetzt. Daher miissen z. B. in Formeln auftretende Abkiirzungen, die in 
Antiqua (d.h. der iblichen Textschrift) gesetzt werden sollen, besonders gekennzeichnet werden (még- 
lichst durch gelbe Unterstreichung). Fettzusetzende Buchstaben sind durch braune Unterstreichung 
zu kennzeichnen. Man vermeide in Formeln und bei mathematischen Symbolen das Unterstreichen 
mit Tinte oder Schreibmaschine; dies wiirde zwangslaufig als zum mathematischen Sinn gehdrig 
interpretiert und daher mitgesetzt werden. Besondere Schwierigkeiten entstehen dadurch, daB 
Schreibmaschinen im allgemeinen keine Unterschiede zwischen 0 (Null) und O (Buchstabe) sowie 
haufig auch zwischen 1 (Eins) und 1 (Buchstabe) kennen. Hier sind unterscheidende Kennzeich- 
nungen unbedingt erforderlich. Die arabische Ziffer 1 schreibe man deutlich als 1 und nicht einfach 
als senkrechten Strich. Das ist besonders verhangnisvoll, wenn die 1 als oberer Index auftritt. In 
englischsprachigen Manuskripten schlieBe man, wenn der kleine lateinische Buchstabe a als mathe- 
matische Bezeichnung vorkommt, die Verwechslung mit dem unbestimmten Artikel aus. 

Durch derartige MiBverstandnisse sind vielfach groBe Korrekturkosten entstanden. In den Kor- 
rektur-Abzigen sollen nur Satzfehler verbessert, jedoch keine inhaltlichen oder stilistischen Ande- 
rungen vorgenommen werden. Nachtragliche (vom Manuskript abweichende) Korrekturen miissen den 
Autoren in Rechnung gestellt werden. 

Samtliche zu einer Arbeit geh6rende (sowohl photographische als auch Kurven und schematische) 
Figuren sind durchzunumerieren. Sie werden getrennt vom Text auf gesonderten Blattern erbeten. 
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